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ABSTRACT 


A one-dimensional numerical model of the urban heat island has 
been devised. It consists mainly in the interaction between a soil layer 
and an atmospheric layer coupled through a heat-balance equation. The 
heat-balance equation includes the various heat fluxes and radiation. 

A constant molecular diffusivity is assumed throughout the soil layer. 
The “second-order" modelling is used in order to close the set of 
atmospheric equations. The equations are transformed into finite- 
difference expressions. The log-linear atmospheric grid is very ap- 
propriate for a one-dimensional model because it gives a good accuracy 
in the low levels with a relatively small number of grid points. The 
initial conditions for the time-dependent model are obtained from the 
steady state. We assess the effect of the roughness height, of the 
geostrophic wind and of the soil thermal conductivity. An increase in 
any of these parameters decreases the amplitude of the diurnal variation 
of the surface temperature. Therefore, these parameters tend to produce 
rural areas warmer than cities. The counterbalancing effect which was 
not included is evaporation which reduces the rural temperature during 


the daytime. 
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PREVIOUS WORK 


1.1 Introduction. 

The observation that the micro-climate of a city is different 
from the surrounding countryside has been made long ago. The most noti- 
ceable difference is observed generally in the diurnal temperature cycle 
with the cities usually being warmer than their outlying areas. This 
is why the climatic influence of the urban area has been termed the "ur- 
ban heet island". The first statistical study of this effect dates back 
to Howard's (1833) study of London's temperature. Since then many other 
features have been associated with the urban effect: changes in surface 
roughness, in soil conductivity, in surface moisture, in albedo, in 
artificial heat input, in air pollution, etc. All these factors inter- 
act in a complex fashion to produce the urban heat island with which can 
be associated not only changes in the diurnal temperature cycle but 
also changes in the precipitation pattern, in fog occurences, in wind 
velocities and in many other properties. Recent observations seem to 
indicate that the maximum temperature difference between the city and 
the rural area happens four to five hours after sunset and not near 
sunrise as could be expected. This feature was observed in Edmonton by 


Hage (1972), and in Vancouver and Montreal by Oke and Maxwell (1975). 


1.2 Statistical Studies, 


Many statistical studies of the urban effect have been made: 
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Duckworth and Sandberg (1954), Landsberg (1956), mitchell (1961), Roden 
(1956) and Bornstein (1968) are a few examples of recent ones. The gene- 
ral results are that the urban area is warmer on the average by up to 2-3 
03 with the maximum difference during the nighttime. In some cases the 
urban daytime temperature was found to be cooler than the rural one. 

The heat island is also generally more pronounced in the summertime with 
a few cities developing their maximim temperature excess during winter, 
This points out that the causes for the urban heat island are varied and 
could in some cases »be minimized or enhanced by local effects due = lakes 
or mountains for example. The wind speed (Sundborg, 1951) and the city 
size (Landsberg, 1956) affect the magnitude of the temperature excess 
with @ smaller value being associated with higher wind speeds and smaller 


cities. 


1.3 Unidimensional Numerical Model. 

The numerical investigation of the urban heat island is recent, 
@eZe, Myrup (1969) was a pioneer in this domain. He used a very simple 
model in which the interaction between a soil layer and an atmospheric 
layer was studied. In Fig. 1, we have a schematic representation of 
Myrup's model. The soil layer contains the minimum number of levelss 
the bottom level is at a distance 2d=0.5m from the surface and is the le- 
vel at which the temperature is held constant; the temperature at the 
level just above the bottom is computed by equating the local change 
in temperature to the diffusion of temperature; the top boundary of 


the soil layer is the earth's surface whose temperature is computed 


from a heat balance equation which contains net radiation, turbulent 


fluxes of sensible and latent heat, and the heat flux from the soil. 
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Schematic representation of the grid used by Myrup (1969). 
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4 
The first atmospheric level is at the roughness height Zo? the next-one 
24 is at the canopy height, and Zo is the top of the boundary layer where 
the meteorological conditions are assumed to be constant. The basic 
assumptions are 
1. Horizontal homogeneity. 
2. Turbulent diffusivities for heat and water vapor are given by the 

turbulent diffusivity for the neareneutral case, 

3. Constant sensible and latent heat fluxes throughout the atmosphere, 
4, Constant relative humidity at the roughness height. | 
5. Temperature, wind and mixing ratio are constant at Zoe 
6. The temperature at a, = temperature at the ground, 
7. The canopy has : unique roughness height. 
8. Constant infrared radiation. 

The most drastic assumptions in his model ares: (a) the lack of 
veeis Glee of the turbulent diffusivities with stability, which is totally 
wrong in the case of stable stratification as shown by the numerical 
results of Mellor and Yamada (1974), (b) the extension of the constante- 
flux layer up to Zo = 300 m, which is a very poor assumption in the stable 
stratification of the nighttime and (c) the assumption that the meteorologi- 
cal variables are constant at Zo» an assumption that Myrup himself found 
quite poor. The most obvious improvement that can be made to his model 
is to add more levels in order to use more realistic assumptions. His 
results showed the maximum excess temperature in the city during the day 
instead of at night as is observed. However, his model, illustrated that 
the following parameters were important: reduction of evaporation in 
the city, increased roughness of the city, thermal properties of the 


buildings and paving materials, and wind speed. Artificial heating 
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did not seem to play an important role in the formation of the urban 
heat island. 

Miller, et_al. (1972) have corrected some mistakes in Myrup's 
model. “Their model is also in one dimension and includes as a new 
feature an albedo model of a city which is being developed by Craig 
and Lowry (1972) and which takes into account radiation trapping between 
the walls of an ideal city. However, this refined model of the albedo 
will not be used in this study because it is not easily applied to a 
real city. , 

By using a more detailed soil layer and more atmospheric levels 
Nappo (1972) has claimed to predict that the maximum temperature difference 
between urban and vive areas occurs a few hours after sunset. However, 
his results should be viewed with caution as he used a constant-flux layer 
of 50 m, which according to Kuo (1968) is not capable of adequately represen- 
ting the temperature variation in the thermal boundary layer. Moreover, 
from his graphs of the diurnal temperature cycle we get the impression 


that the temperature variation would be different from one day to the othor. 


1.4 Two-Dimensional Models. 

Some recent studies have extended the model to a second di- 
mension, such that horizontal advection is possible. For example, 
McElroy (1972) has studied the effect of land-use strategies by alter- 
nating rough and smooth surfaces in his model, with the purpose of trying 
to derive benefit from the urban heat island. However, he used a constante 
flux layer of 15 to 35 m, which is certainly too high for the actual overe 


night constant-flux layer. 


Are there many advantages in going into a two-dimensional 
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model? Part of the answer lies in the way in which the two-dimensional 
model treats the discontinuity in roughness height between two adjacent 
horizontal grid points. In general the models assume a fixed roughness 
height at each grid point. Therefore, at each grid point we assume ine 
finite homogeneity in the transverse direction and spatial homogeneity 
of the order of the grid-point spacing along the grid-point orientation 
direction. From studies of the constant-flux layer, it is clear that 
discontinuity in roughness heights causes the formation of an internal 
boundary layer as illustrated in Fig. 2. Above the internal boundary 
layer the turbulence is characterized by the upstream roughness height 
and not by the roughness height below it. The ratio of change in rough- 
ness height to fetch is approximately 1/20 on the basis of the wind 
profile, according to Angle (1973). Im many two-dimensional models 

the horizontal grid-point spacing is about 500 m, and generally the 
change in roughness height is at most 1 m between each grid point. If 
we assume that the change in roughness occurs abruptly in the middle of 
two grid points, we get an internal boundary-layer height of 12.5m 
over the second grid point as shown in Fig. 3. This means that all 
levels above 12.5 m are almost undisturbed yet some authors have assumed 
a constant-flux layer of up to 50 m, which means that they have assumed 
the wrong value for the lower boundary condition and that they have in- 
troduced too soon the influence of the change in roughness height at 
higher levels. Therefore, the two-dimensional models might not be able 
to represent accurately the transition from smooth rural area to rough 
urban area. However, the two-dimensional model has more versatility 


certainly because it can include vertical winds and divergence. 
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transition layer 


Figure 2. Schematic illustration of the horizontal variation of the 
depth of the internal boundary layer (h) and of the height of the tran~ 
sition layer (h'). The roughness height changes abruptly from a to b 
at Xo The region above the transition layer is influenced only by 


Z.~ =a, whereas the region below the transition layer is influenced only 


ro) 
by 2,= be In the transition layer both a and b are influencing the 


flow. 


xy discontinuity Xo x 
me Peg) 
Figure 3. Diagram of the horizontal variation of the depth of the in- 


ternal boundary layer for a two-dimensional model in which the discon- 


tinuity in roughness height happens between two grid points. 
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1.5 HighoreOrder Closure Scheme. 
All the models discussed so far have used K-etheory in which 
the stress components are related to the gradient of the mean wind 


components by, 


aw e- x, 24 Gis 1) 
w=-K, ou (1.2) 


where, Kn = eddy diffusivity coefficient for momentum 
U and V = x and y components of the mean wind (horizontal) 
uw and vw = x and y components of the friction velocity u, 
related to the stress T and air density P by: 
uy, = (t/p)? (1.3) 
Arbitrary assumptions have to be made in order to define Kine 
Close to the ground in near-neutral stability K, has been observed 


empirically to behave like, 


xv 
Ky ~ key (z + 2/3 | (1.4) 
where, ky = 0.35 = Von Karman's constant ( Businger, et ale, 1971) 
Ds roughness height 


z = height above the roughness height. 

K-theory has been extended outside the constant-flux layer by 
assuming an empirical variation of K, with height. The near~neutral 
stability condition can be removed by multiplying K,, by an empirical 
“stability function". All of thase extensions do not appear to be ona 
solid theoretical basis and depend rather strongly on the experimental 
data used to derive them. Therefore, while K-theory has the great ad- 


vantage of simplicity, it can be very poor on occasion; for example 
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a countergradient flux has been observed during some experiments (Eski- 
nazi and Zrian, 1969) in locally-embedded regions where the stress and 
velocity gradients are small and this cannot be explained physically by 
a gradient theory like K-theory. An alternate approach in numerical 
modelling of the boundary-layer is the "higher-order closure" scheme, 
being developed, for example, by Donaldson (1973), Mellor and Yamada 
(1974), Meroney (1974) and many others, Mellor's and Yamada's approach 
seems the most promising for the study of the boundary layer. They have 
ordered various levels of models in terms of their simplifying assumptions, 
and they have imposed boundary conditions near the surface in accordance 
with experimental results. The weakest point of these theories is that 
the constants used in modelling various terms in the equations are ill- 


defined due to the lack of precise observations. 
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CHAPTER II 


DEVELOPMENT OF THE ATMOSPHERIC HQUATIONS 


2e1 Basic Equations. 


2e1ei1 Perfect Gas law. 

The equation of state that will be used is the perfect gas law 
which is accurate enough in comparison with the uncertainties associated 
with many turbulent variables. However, as this model contains water 
vapor in variable amounts, we will use the virtual temperature concept 


in order to keep R, the dry atmosphere gas constant, in the equations, 


te LNG 
PORT, (2.1) 
where ™~ indicates the instantaneous value of a variable 


? is the pressure 


ie is the virtual temperature related to the absolute 


“~/ 
temperature T by, 


“vv 


T, =T (1 + 0.6078 4) (2.2) 


wheres q is the mixing ratio for water vapor. 


2e1e2 Continuity Equation. 
The continuity equation represents conservation of mass in 
a compressible atmosphere. 
soca 4) = 0 (2.3) 
where , is the wind component in the j direction, and 
x is the distance in the j direction. 
10 
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21.3 Momentum Equation. 
We will use the Navier-Stokes’ equations to relate the stress 


~~ 


Cre and the rate of strain q.. These equations require the following 
assumptions: 
1) There is a linear relation between stress and rate of strain. 
Ts 5 ==) Si +h 
“nw IU 
o- 
where 43 = KG sui SS 3x, 
and Ks jk is a fourth order tensor whose properties are defined 
by the next two hypothesas. 
2) The molecular structure of the fluid is isotropic. 
3) The divergence makes no contribution to normal stresses. 
4) The variations in the first and second coefficients of viscosity, 
respectively and Kee are neglected, 
If we add the Coriolis force and gravity the equations of moe 


tion become: 


fy Puy erie ee Dit aes re 
0st Sm Cee “ax “Ox, Tig 7 P S7PEk ja £5 Ue (204) 
where the stress equation is 


and where: i, j, k are indices of value i to 3. The repetition of an 
index inside a term implies summation over all the values, 
+1 if i, j and k are in cyclic order 
sk = 4=1 if i, j and k are in nonecyclic order 
O if any of the indices are equal 
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201.4 Energy Equation. 

Following Godson (1958), we start with the equation for the 
entropy of the dry air - 

dS, =C, dln T~- Rd In (p-8) (2.6) 

where @ is the partial pressure of water vapor and 

Ch is the heat capacity at constant pressure of the dry air. 

A similar equation is valid for S, the entropy of the 


water vapor: 


dS, = C dent = Rody ine Caw 
where Cp, is the heat capacity at constant pressure of the water 


vapor and R,, is the gas constant for water vapor. 
The following relationships are useful (Godson, 1958): 
€=R/R, 
Cp, = 8 Cp/(7E) 
6 = p a/(qt€) 
p- 6 =ép/(qté) 

We combine (2.6) and (2.7) and form the equation for Ss, 
the entropy of the mixture of 1 gram of dry air with q gram of water 
vapor. If we neglect some small terms involving (aq), we obtain 

(1+9)a5=(c, +40, ) din T - (R +qR,) d In 9(2.8) 
We use the definition of virtual temperature and the perfect 


gas law, neglect some small terms including the variation of the heat 


capacity due to water vapor, and obtain after some manipulations 
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we obtain an equation stating that the change in entropy of the dry 


air is due to molecular thermal diffusion, 


aT eS de (2.10) 


A similar equation holds for the change in entropy of water 


vapor due to molecular thermal diffusion, 


Gace 
(ol < EASON OL) 
pyt 48 pcveaioxt (2,14) 


Again we compute the change in entropy for the mixture of 


one gram of dry air with q gram of water vapor, 


pac tad Gi 2 ro GV Was ee Wome 
Gu Se Tt + Ty ore ee ¥) (2512) 
Pty at Cs OX, 9X; sal Eas) Srey © 


In the last equation we have neglected two small terms of 
the order of (AG)* and (AT AG) with respect to AT and T Aq for the 
other terms. Using the equation for water vapor conservation (2.14), 


we combine (2.12) with (2.9) in order to get the energy equation, 


ov 


dty SORE 90 eae ey 
Cy Cos + Belo Gy Ty)= 2 + Gy SP + Rete SED (2.13) 


201.5 Water Vapor Conservation Equation. 


We make the hypothesis that the change in the mixing ratio 
following a parcel of air is due only to molecular diffusion; therefore 


we neglect condensation and eveee Gums 


“~~ 


aN Sas ea 
pres Rath (2dr Gok 2y.282)7 oe CW 28 ox ean 
J J ra 
In this equation we have implicitly assumed that the physical 
processes responsible for the diffusion of momentum are the same as 
those that cause molecular diffusion of temperature and of water 


vapor. Consequently, we can write 


~s ine 
vi = k/Cp 
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2e2 Assumption of a Nearly-Adiabatic Atmosphere. 
If we consider the atmosphere as being in a state slightly 


removed from the dry adiabatic atmosphere at rest, we can express the 


variable A as being the sum of its adiabatic value A, and a deviation 


A’, 
Pp =p, +P! (2.15) 
Gj= 0 + uy, (2.16) 
ew ae Ne (2.17) 
P= C3 cas (2.18) 
T= Ty,+ Ty (2019) 
q =0+q! (2.20) 


2e2ei1 Perfect Gas Law. 

Substituting (2.15) and (2.18) in the perfect gas law (2.1), 
we get after separating the equation into an equation for the adiabatic 
state and an equation for the departures from the adiabatic states 

te = Po R Ty, (2.21) 


p’ = RCP, TL +P" Tyt P* 4) (2.22) 


Perec Continuity Equation. 


Similarly, the continuity equation (2.3) becomes, when 
combined with (2.18) and (2.16): 


en | 
: = =. (2.23) 
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20203 Momentum Equation. 


We now make the additional assumptions that 
2 ‘ 
ou ASS Po MK fo Ay Sc a, (2.25) 
and substitute (2.15) to (2.19) into the momentum equation (2.4) 


making use of (2.22) and (2.25) in order to simplify the terms. 


9 Po _ 

Sz, 8 2426) 
> : su oe } > yu! yu 
Pea te PS GO. Fe) + 
J L J 5 3 
Ae tye Un | D2) 

LJ ayes S a (2227) 


20204 Energy Equation. 


We now make the substitution into the energy equation, in 
which the thermal diffusivity has been rewritten in terns of M + The 
following additional assumption is needed in order to simplify the 


equation: the heating of the fluid due to motion ( ou + us 57] = 


O (Mach)*) is neglected. 
The entropy equation (2.8) is applied to the basic isentro- 
pic atmosphere and becomes, 
a tye 2 So eh CARGoncDSt” G¥ac Gite (2.28) 
If we take the local time derivative of (2.28) and note that 


the entropy of the basic atmosphere is time independent, we haves 


T p 
oe es CD oo v. a (2.29) 


The equation for the energy of the basic atmosphere as de- 


termined by the substitution of (2.15) to (2.19) into (2,13) is, 
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a i 
ie Oe YPo Vo _ 
s—)) = UGE y Seem nmares he (2.30) 
The last equality results from (2.23) and (2.29). Therefore, 


we are left only with the energy associated with the diabatic part of 


the atmosphere, 


T 
SE ee ey eet ano 
are Pov st ‘g oy ox ) © OX. O x5 +} dx; ) Cy (2531) 


20205 Mixing-Ratio Equation. 


When we introduce the adiabatic reference state into (2.14) 


we have, 


Os SemitaP obsess 


oot 


dq’ _ 3 T dq? " 
x; ied eter (2.32) 


202.0 Simplified Continuity Equation. 


From the perfect gas law equations (2.21) and (2.22), we 


can form the ratio of the adiabatic deviation of the pressure to the 


pressure of the adiabatic atmosphere, 


’ ’ ’ 
peso piesa e (2533) 
Po Ty, Po Po Vo 


and sel Oe » the equation reduces to, 


Asp! ¢< Po 
es 

pit civ aie! (2.3!) 

p 


The changes in p, are of the order of ( p, ut) so that, 


aa Oph) =0Cy macn) = Mel ea, 


Ty Po 
where ¥ is the specific heat ratio of the atmosphere = Co o/ Cv, and 
is the heat capacity at constant volume of the atmosphere. 


et =" Ty / To (2.35) 
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17 
Combining (2.35) and (2.24) and neglecting second-order 


terms, we obtain 


GeBjedaret peo tuk — pe he 


aX, Po eX PS i Ho § (2.56) 


For the atmospheric processes in the boundary layer this 
equation can generally be simplified to, 


ou 
> X4 = (237) 


2e2e7 Resulting Set of Equations. 
If we replace Ty by its approximate expression, 
besiege (14+ 0.6078 q?) (2.38) 
we can then easily derive the equation for T* which has the same form 
as the cnes for ee and q’. Therefore, the atmospheric equations in the 


quasi-adiabatic atmospheric model ares 


Po via Ge R Tyo (2.39) 
Sania (2.40) 
oY a itl ao bile ton eee ne ee 
Osean tot J Serre ody ly meer tHe Szioa; 4 
a f wy. (2.41) 
oT? 1 
Po It +P, uf i is 3X5 ae (2.42) 
oq' ode Lye 
Poseatipe 34 " “Wo = x, q (2.43) 
Tue 
Qo _ ta (24) 
Su 


=0 (2.45) 
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Z2e3 Ensemble Averaging. 
The next step is to separate each variable into a mean value 
and a departure from the mean. We will use the convention that a vae 
riable written in uppercase means an averaged quantity whereas if it is 


in lowercase it represents a fluctuating quantity. 


we =U, + uy (2.46) 

tT? =Q+6 (2.47) 

tT =@) +0, (2.48) 

q?ye each (2.49) 

elt (2.50) 

Mis A (2.51) 

V =)A/0 = kinematic coefficient of viscosity (2.52) 
= Pp (2.53) 


2e301 Equations for the Mean Quantities, 
The equations for the mean quantities are obtained after 
inserting (2.46) to (2.51) into the equations of Section 2.2.7, and 


then taking ensemble averages of the new equations. 


Quy 
Ey = 0 (2.54) 
ih ae _ 1 2p, 84D 
rt sepa be Beij) aaeahat ice ah 6 Aeoge Stet 
rt U (2.55) 
axed 
> 5 = ene 
ty "xe ( u.© + U, ) By (2.56) 
Yq 0% a 
yt eT ( U, Qt Uy. wm qd y) =¢4 ——= 5 x2 (2 57) 
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2e3e2 Equations for the Turbulent Quantities. 


The equations for the turbulent quantities are derived from 
(2.41) to (2.45) in which the decomposition into mean and fluctuating 
parts has been made. From these equations are subtracted the appropri- 


ate equations for the mean quantities as given in the previous section. 


ue 
=! =0 (2.58) 
J 
dus vests 
Sor tse ie ty tM US to Wy = BER) +Ega fue 
BAe Sie ys (2.59) 
Co oX5 abe dxe 
2 
oes era Oty t+uG)- yE acess (2.60) 
Ak 
2 
Cicer: (Ua) 86 a tae Ue g) Be a (2.61) 
5 k 


20303 Equations for the Second-Order Moments. 


The equations for the mean quantities involve the second- 
order moments for which equations will be derived from the basic set of 
(2.59) to (2.61). The equation for a, Uy is obtained by changing the 
subscript j to i in equation (2.59), by multiplying the resulting equa- 
tion by Vy and adding the new equation to (2.59) which has been multi- 
plied by u,. The final equation for Uy U5 is then obtained by taking 


the ensemble average of the last equation. 
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Mek ‘gpk otis 
- (8; a Oo SF Bs U5 Us Oy Me raya am, ea lees 2 ace (2.62) 


Similarly the equation for a5 is obtained by multiplying 
(2.60) by u, and (2.59) by ©, adding the two resulting equations and 
finally taking the ensemble average. Because the equations for the 
other moments are found in a similar fashion, we will simply write 


down all the other equations without any further explanation. 
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2+ Modelled Equations. 
Equations (2.62) to (2.67) are the equations for all the 


second moments needed to solve the equations for the mean quantities, 
Unfortunately, these equations involve higher-order moments as well 
as covariances with pressure and many other covariances. In order to 
close the system we must find equations for these new covariances. 
One method would be to derive formally equations for these variables, 
but that method would introduce new covariances for which equations 
would be needed. This course of action does not look very promising 
because we would obtain an infinite set of equations always with 
fewer equations than variables. An alternative method consists of 
modelling each of the new covariances in terms of know variables 
for which an equation exists. We will follow Mellor and Yamada's 
(1974) modelling assumptions, There are some general guidelines 
that can be followed in that type of modelling: 
1- The tensor properties of the term to be modelled should be 
preserved if possible. 
2= The dimensions should always be respected. 
3- The sign of each modelled term should be selected in such a 
way as to avoid negative scaling lengths or velocities. 
4- We should try to use the simplest model possible which obeys 


guidelines 1-3 except if results indicate that a more complex 


model is needed» 
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24.1 Energy Redistribution Term. 


The name “energy redistribution term” was used by Rotta 
(1951) when he noticed that this term partitions energy among the 
three energy components while not contributing to the total. Monin 
and Yaglom (1970) discuss Rotta's work and state that the physical 
justification for the modelling of this term lies in the Poisson's 
equation for pressure which is obtained after taking the divergence 


of the Navier-Stokes’ equation without external forces, 


Apt = (A x) (2.68) 


4) ares 

If we apply Reynold's decomposition for ensemble averaging 
to (2.68) and subtract from it the equation for the mean pressure ob- 
tained by taking the ensemble average of (2.68), we are left with the 


equation for the fluctuating PESeoure 


>2 (Us Weer Us uy tu, us -U;, Uy ah (2.69) 


ine See ete 
This equation applies to a point element of a fluid, whereas 
the pressure fluctuation in our equations is representative of a verti- 
cal column and p therefore would result from integration of Ap over a 
vertical colum. We will not carry out this integration here because 
we are interested only in the form of the terms that can be extracted 


by the inclusion of (2.69) in the energy redistribution term. 
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3 U: 
when we average this product we get terms of the form oe 


where Cy jkm is a fourth.order tensor yet to be determined. The other 


terms will involve products like bP Uy vy Quad which are likely to 
dX 9X4 9X4 


give rise to terms in Uy Uy. Therefore the energy redistribution term 


can be modelled by an exuression like, 


Us gu OU 
6, 55 ‘3%, hace Cs jum "km tC ijk are (2.70) 


Rotta (1951) then assumes that the constitutive coefficients 
are isotropic. That assumption has the advantage that it respects 
the tensor propertiss of the energy redistribution term and gives a 
relatively simple final answer, 
C5 sim = %4 S25 Siem + 2 ste Spm + 3 Ste Sam (2.71) 
Caen = Oh S44 Sim +53 See Fam * SF je Sam (2.72) 
Substituting the isotropic tensors (2.71) and (2.72) into the 


terms on the right-hand side of (2.70), we find 


c OU pee Riel 
and Chm Se = C4 > = + os 2d Saye (2.74) 


The equation of continuity has been used in (2.74). Due to 
the symmetry of the terms involving C4 and ss) it seems justifiable to 
postulate that these coefficients are equal and, therefore, thats 


9 U,, yu U 

. i, 0; 
ijkm soo 

If we set i=j (1=1, 2, 3) in (2.70) and add the three resul- 


) (2.75) 


ting equations we get, 
Seno 2 2 2 JU 
cess ee + (Cote +2C = 2.76 
pean +? 3 C, @ (Co 3) 8 2 ix, (2.76) 
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2= Uz, Us is proportional to the turbulent kinetic energy. 


where 9 


Using the continuity equation, (2.76) reduces tos 


Ca reC 
Cy ad = (2577) 


Now we compare the dimensions, 


Es See + eee ab -[F)-[r ¢] (2.78) 


lie, “| as (c, + C3) A (2.79) 
[ios ¢ a + al ))}- \c3 z | | (2.80) 


where L, HM, T, and U indicate quantities with dimensions of length, 
mass, time and velocity, respectively. Therefore, on dimensional 


grounds, 


(Cp +C4) = 4-3 (2.81) 


Le) 
AS 
— 


CC aca (2,82) 
where sir is a sealing length that has to be determined empirically, 
and C a constant that Mellor and Yamada have determined to be 0.056. 
Now we have the sign to consider. The easiest way of determining the 
Sign is to examine the variation in the vertical because the vertical 
gradients are normally the largest. We will first determine the sign 
of the energy redistribution term. Let us imagine a horizontal velo- 
city perturbation u with a vertical profile as showm in Fig. 4 extene 
ding over the range 2 h centered at 2,- From the continuity equation 
we would expect the vertical wind to compensate for the increase in u. 
Furthermore, we imagine the ideal situation in which there is no mixing 
of the air above and below an and in which half of the compensation 


comes from below and half from above. From Fig. 4 we can infer that 
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Figure 4. Shematic representation of the vertical distribution of the 
variables in the energy redistribution term. The zero value is indie 
cated by a solid line parallel to the 2 axis, and the dotted lines 


represent the behavior of the variables. 
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air is removed around Zt h and transported to levels near Zoe There= 
fore, there is a decrease in mass between want h/2 and z,+ h creating a 
negative pressure disturbance, and there is accumulation of mass below 
that layer down to Z, Which causes the pressure disturbance to become 
less negative. As there is no mass transfer across Ze We expect no 
pressure disturbance at Ze if we neglect the dynamic effect caused by 
the downward vertical wind. Similar reasoning applied to the region 
below Ze gives upward vertical motion and a positive pressure disture- 
bance. By comparing the u and p vertical profiles we see that the proe- 
duct p A is always positive. This idealized situation clearly shows 
that we expect the energy redistribution term to be positive when we 
consider the vertical variation of the horizontal wind. 

If we use (2,82), the first two modelling terms in (2.76) 
becone 377 (uj uy - e“/3). As we know e* is proportional to the tur- 
bulent kinetic energy and is always positive. uwis generally the 
most important covariance included in wu, Us and can be related to the 


J 


gradient of the mean wind through a diffusivity equation, 


K,, is generally positive and the gradient of the mean wind is 
always positive in the nearly constanteflux layer. Therefore, the first 
two modelling terms are always negative which implies that we need a ne- 
gative sign in front of them, The modelling term C = is always posie 
tive so that its sign agrees with the sign of the term that is being 


modelled, Grouping all our findings, the final result is, 
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This term has been modelled by many boundary-layer workers 
with exactly the same terms but with differences in the values of the 
coefficients. The coefficient C has been assigned a value of 0 by 
Donaldson (1973) and a value of 0.4 by Deardorff (1973). These can be 
compared with a value of 0.056 chosen by Mellor and Yamada. These dif- 
ferent values show that the determination of the coefficients is not 
unique but depends on the experimental results used to derive these as 
well as on the values assigned to the coefficients of the other model- 
led terms, There is a closer agreement on the value of the length 
scale 3/4. In the constant-stress layer Donaldson has chosen 2 value 
of 0.7 2 which compares well with the value of 0.822 selected by Mellor 
and Yamada, Deardorff's results are not directly comparable because 
his model has to be applied outside the constant=stress layer. Instead 
of using a mixing-length approach he uses a Scaling length proportio=- 
nal to the grid spacing which he thinks is an appropriate length for — 


modelling sub-grid seale phenomena. 
ply oe 
ea ox, 


This term is similar to the energy redistribution term and 


204.2 Pressure-Temperature Correlation Term: 


by similar arguments we can hope to model it by, 


Pe, Sere ae _ (2.84) 
Po 2X5 j OX, 


We match the dimensions, 
0¢| . |U 
ee =A [Evel (2.85) 
[p a6 |= [p vel (2.86) 
= E ve | (2.87) 


Mellor and Yamada have put D*=0 and used f , as the appro- 
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priate scaling length, 
e 


” 3R2 
The magnitude of the modelled term is now defined and its sign 


D (2.88) 

- will be determined in the discussion which follows. Again we will 
examine only the vertical variation because the gradients of most of 
the atmospheric variables are almost vertical. In thse case of the ener-= 
gy redistribution term we used the fact that the gradient of the mean 
horizontal wind was always positive in the lower boundary layer, 
However, this simple situation is changed in the case of the vertical 
gradient of the mean temperature. Here a positive gradient represents a 
Stable situation whereas a negative gradient is an unstable situation. 
In a complete discussion we should analyse both situations, but here 

we al be satisfied to find the sign of the modelled term in the 
unstable case, Fig. 5 represents an unstable temperature profile as a 
function of height. We will use the method of the air parcel familiar 
to all meteorologists in order to get some physical insight into that 
term. The justification for using that method is that we are looking 
at small elements of the atmosphere which will conserve their identity 
during a short period of time. The dissipation term is very ef= 
ficient for small atmospheric scales and therefore this method is valid 
only approximately and for very short distances. If an air parcel is 
lifted adiabatically and without mixing with the environment from 24 

to z> the air parcel is warmer than the environment by O(2,)- Mz.) 
=@,. If the air parcel had been lifted from z, to 23 the tempera- 
ture difference becomes ()(24)- W(z5) = 0,7 94. Therefore an unstable 


situation is very likely to give rise to 2250, In the atmosphere, a 
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Figure 5. Illustration of the air parcel method. The solid curve re# 
presents the mean potential temperature of the environment for an 


unstable situation. 
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warmer air parcel implies a lower pressure and, therefore, this warming 
is associated with a negative pressure disturbance. The result in the 


case of an unstable lapse rate is that, 
29 
Dig (2.89) 
The modelled term in an unstable stratification is, 
we =- Ky ABy0 (2.90) 


K, is always positive in the lower boundary layer when the 
lapse rate is unstable and, therefore, we need a negative sign in front 


of the D term. The final result is, 


a é 


pa Seca Bee ae (2.91) 
Donaldson (1973) has kept the same term with the scaling 

length &, = 0.7 z in the lower boundary layer; this compares very 

well with the value of 0.82 z that Mellor and Yamada (1974) have 


used. Deardorff (1973) added one more tern, 


(2.92) 


This term is valid only for isotropic turbulence and corres~ 
ponds to a gravitational compensation effect. Deardorff (1973) admits 
that he does not know how rapidly this effect becomes invalid as depar- 


tures from isotropy become more pronounced. 


2.4.3 Pressure-Water Vapor Correlation Term: aa 
cr) 


Again we postulate that we can model t. = term by an expres- 


sion like, 
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mae =EBUya (2.93) 
fe) 
The dimensions are, 
| ay es E UQ (2.94) 
Paes L ] 
[Bua | = [gua] (2.95) 


Therefore, 


ae Ane) 
ree em (2.96) 

where 1 is a scaling length to be determined. The sign is derived by 

a demonstration analoguous to the one for temperature. When we consider 

the case of upward water vapor flux the arguments are exactly the same 

as in the preceding section if we replace the petit al temperature 

by go =Q+ 0.6078 q and consider the changes in Ge to be caused only 


by the changes in q. Therefore, a negative sign is needed, 


Pods. 2, 
Ee 315 U5 4 (2.97) 


204.4 Dissipation Terms. 


The dissipation terms are the ones involving the viscosity 
coefficient and the product of the gradient of two turbulent quantities. 
Following Kolmogoroff (1941) we hypothesize local small-scale isotropy. 
Therefore, all the dissipation terms will be modelled in terms of the 
simplest possible isotropic covariances. 


algenhess AY 
sine Fe oe. (2.98) 


The dimensions on the left side are aS, therefore LF) = 
[u/u). A will be used as the length scale for the dissipation terms, 


with e as the scaling velocity. The positive sign is needed if we 
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want destruction of covariance and not production of it. The dissi- 
pation of covariance involving any two different wind components is 
zero because it would require an anisotropic modelling term of the form 
Wy ye The only isotropic covariance involving the wind is ef = Wy, Wy 
as was used in (2.98). Therefore it seems reasonable to model the dise 


sipation term by an expression like, 


yet 25g 2 2 0? 


— == 209 
39%, 3 Ay O83 ea 
The other terms are modelled similarly, 

eee Ae 

242928 =26 9° (2.100) 
dX, 9X > 
seeeecn Re 

ood 24 = 26 (2.101) 
IX, 3X, \3 
CR ee Oe 

2¢ 20 =2e (2.102) 
> XO Xp Ay 


The remaining two terms would require modelling terms like 
aye and us which are not isotropic. For this reason they are put 


equal to zero. 


Cah ea talon), Py (2.103) 


24.5 Diffusional Terms. 

Most workers use a flux-gradient model for the various 
triple-moments terms. This implies a turbulent redistribution that 
leads to smoothing and hence the terms are often called "diffusional", 


- The formula proposed by Mellor and Yamada (1974) is, 
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That formula appears appropriate and we determine that LG] = 
[L U] from a consideration of the dimensions. We choose A as the sca= 
ling length for the diffusional terms, keep e as the scaling velocity, 
and give a negative sign to the expression in order to be consistent 
with the gradient form in which it is written. We have finally, 


dh sleek gu, Us 
U, Uz Us = = ed aes at k 2 (2.105) 


More complex expressions have been derived from a more phy= 
sical basis. For example, Hanjalic and Launder (1972) have modelled 
the pressure correlation terms by terms similar to the ones derived in 
sections (2.4.1) to (2.4.3) and then have derived the equations for the 
third moments. They have applied the quasi-enormal assumption to the 
fourth moments and have neglected the remaining terms. The final re- 


sult is a term like, 
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U,~ Uy Uy = G! (au, Uy 
(2.106) 
Wyngaard (1973) has shown that the assumption of free con- 
vection predicted that both w andd w*/} 2 were positive. These pre- 
dictions are partly supported by the results of the Kansas experiment. 


Therefore an expression like, 
weg (se = (2.107) 


would require a positiva value for G and not a negative value as 
in our model. Therefore, the modelling of the triple correlation terms 
will prove to be very poor in all cases where free convection dominates. 


However, we are modelling forced convection and Wyngaard's finding 


might not apply there. 
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Using the same assumptions the remaining triple covariance 


terms are modelled similarly, 


2 
u,0% == eA, (oe (2.108) 
ee 2% 6  atye 
Wy. 450 7 CGE Near rp ea x, (2.109) 


ae ery yx; Sas ee (2.140) 
2 »@ 

Uy oS eA. DX, (2.111) 

Uy, a9 =o ods = 4 (2.112) 


Again we expect these expressions to be valid in forced 
convection as is the case in the surface layer with relatively strong 
winds. In order to model correctly these terms in the freeeconvec- 
tion limit we would need supplementary terms probably involving 


buoyancy. 


2.4.6 Pressure Diffusional Terms. 


Hanjalic and Launder (1972) have found experimentally that 


these terms were small. Therefore, we shall neglect these terms, 


Ma po = Pde 0 (2.113) 
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Some authors have proposed models for these terms. For 
example, Donaldson (1973) has used a model in which the terms are 


nonproductive as in, 
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where \’ is a scaling length. However, Donaldson (1973) reports that 
Taulbee (1973) had to use a small negative value for A’ in order to 
obtain better agreement between computations and experimental results. 
This means that it is better probably to neglect the pressure diffusion. 
al terms until we have sufficient experimental results to suggest 


a physical basis for the modelling of these terms. 
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2e5 Modelled Equations. 
2e5e1 Modelled Equations for the Second-Order Moments. 


We insert all the modelling terms into the second-order 
moment equations (2.62) to (2.67) and neglect the molecular diffusive 
terms as well as those involving the Coriolis force which are a few 
orders of magnitude smaller than the other terms. The modelled equa 


tions for the second-order moments ares 
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The Lagrangian derivative is defined as, 
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2e5e2 Modelled Equations for the Mean Quantities. 


The equations for the mean quantities are not directly affec- 
ted by the modelling and, therefore, we will use equations (2.54) to 
(2.57) in which we neglect the small molecular diffusion terms and 


replace the Eulerian time derivative by the Lagrangian one. 
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2.6 Length Seale Specification. 
Up until now the scaling lengths used in modelling J, , me \, 


and the constant Cy have not been specified. If we accept the results 


obtained by Mellor and Yamada (1974) we get the following relationships, 


p, =J,.=45 = 0.78 (2.12!) 

Oy =), aA =) = =A = 0,23 (2.125) 
\, = 15.0 (2.126) 

N, =1, =A, = 8,0 (25127) 


We have assumed that water vapor behaves like temperature with 
respect to turbulent properties, so that their respective length scales 
are identical. f is the fundamental length to be determined. We will 


use Blackadar's (1962) formula, 


(z+z,) 


6128 
Prescn (2.128) 


ne 4s a constant which fixes the size of the largest eddies 
and whose value depends on the stability of the atmosphere. That fore 
mula holds for the two known asymptotes of the length scales close 
to the ground it behaves like jz k,(2tzo)» and much above Zo 2h 
reaches its maximum value ip We must specify i and to do SO, we 


adopt the arbitrary formulation of Mellor and Yamada (1974) which is, 
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(2.129) 


Equation (2.129) gives a reasonable number of the order of 
30 m in the case of neutral stability. Di increases with stability 
and decreases with instability. Finally, C, is given the value 0.056 


following Mellor and Yamada (197%). 
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2.7 Ordering of Terms. 


Mellor and Yamada (1974) have described a method of ordering 
terms which works quite well. The method consists in separating the 
isotropic part from the anisotropic part in all of the equations. That 
task is simple as only equation (2.114) consists of a mixture of the two. 
Equations (2.115) to (2.116) contain only anisotropic terms and equations 
(2.117)-(2.119) contain only isotropic terms. Going back to equation (2.114) 
we note that the isotropic part is the energy equation obtained by the 
addition of the three equations in which i=j, for i1=1,2,3. Therefore, 


the equation for the isotropic part is, 
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The anisotropic part is found by subtracting j, /3 times 


equation (2.130) from (2.114). 
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The purpose of that exercise is to involve the ratio 2/N 
in the equations. From neutral experimental data Mellor and Yamada 
have shown that this ratio is generally of the order of 0.05 to 0.10. 
This gives an objective method for ordering the terms by neglecting the 
terms which are smaller by a factor of the order of M/A. Table 1 is 
the result of such orderings; the first lines in each block are the equa- 
tions for the second-order moments in which (2.132) to (2,135) have been 


used; the second lines in each block are the order of magnitude of each 
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Table 1. ori of Terms, This table is similar to table 1 in Mellor 
and Yamada (1974) page 1794. 

- The first lines in each block are the equations for the second-order 
moments in which the departures from isotropy a, b and c are used. 

- The second lines are the order of magnitude of each of the terms. 


- The third lines are the order of magnitude of the dominant terms. 
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lity 
term as given by equations (2.136) to (2.146); finally the third lines 
in each block are the order of magnitude of the dominant terms in each 
equations these lines are obtained when we make use of the relationships 
between the various terms, which will be derived in the following dis- 
cussion. 

As a first step, we will discuss the simple case of neutral 
flow by neglecting buoyancy which is done by giving a value of zero 
to ge We assume that the production and dissipation terms are dominant 
in all of the equations (2.147) to (2.153). The production terms are 
the ones involving the gradient of the mean quantities, whereas the 
dissipation terms are the ones containing e times the second-order 


moment considered. We apply this assumption to (2.147) and (2.148), 


ae’ U, = ee (2.154) 
oe U =a off (2.155) 
We combine the last two expressions to get, 
a? = S/N (2.156) 
vu, =4/(ah) (2.157) 


The domination of the production and dissipation terms implies 
that the right-hand sides of (2.149) and (2.150) are the most important 
terms, therefore 

e 2 b0,= o9*/N (2.158) 
eVo Qe oV/N (2.159) 
Using this assumption in the other equations we have, 
o°@, =e 9 v/f (2.159) 
e% Q =o Vo/f (2.160) 


The last four equations yield the following, 
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be = Ayn (2.161) 
cf = hin (2.162) 
@, = P/(bN) (2.163) 
Q, =W/(e A) (2.164) 


Obviously we have a=b=c. If we allow buoyant energy produc~ 


to be of the same order of magnitude as the dominant parameters, we 


get, 
2 
eg? = eo ® 
be bA (24165) 
ce) 

a 

PMc 
0.61 Ae c IN (2.166) 


The ordering of the terms has not been completed yet because 
the diffusion and advection terms involve a wind scale U and a length 
L whieh are undefined. Mellor and Yamada (1974) have made a 4mlevel 
classification. The level 4 model is the most complex and keeps all the 
terms. The other levels are obtained by assigning a value to the ree 
lative order of magnitude of the diffusion and advection terms. For 
example level 1 is the simplest model in which the diffusion and ad- 
vection terms are neglected as well as other terms of order aaC LA)? 


In the next section we will discuss in detail how level 3 is obtained. 
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2.8 Level 3 Model. 


We consider the equation for the kinetic turbulent energy oe” 
and look at the dominant component of the diffusion term which is the 
vertical component. We assume that the vertical derivative is of the 


order of 1/ A so that the diffusion term becomes, 
Vr ay 282 EO AA. ys Bee! 
ST Ones oe 1G CC nied) © OF aD (2.167) 


where the last equality comes fron WA = @ (a%) which seems to be 
justified from experimental data and numerical boundary-layer model- 
ling ( Mellor and Herring; 1968,1973). Level 2 keeps (2.167) intact. 
However, this relationship is quite uncertain and it is possible that 
the diffusion terms without being dominant, are more important than 
(2.167) implies. Therefore, in level 3 we replace a“ by a in (2,167) 


and get for the order of magnitude of the diffusion term, 


2 
oe =a ce (2.168) 


Accordingly, the ratio U/L is defined and can be used in 
(2.147) to (2.153) to get the final estimate of the order of magnitude 
of each of the terms. Next, we multiply every ordering term in (2.148) 
by a, in (2.151) by b and in (2.152) by c and neglect all the terms 


of O (a2), G(b@) ana O (c“). The final result is, 
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2.9 The Boundary-Layer Approximation. 
We will restrict ourselves to the boundary layer in which 
the two following approximations can be made. 
1) The vertical gradients are very much larger than the horizontal 
gradients, 


ioe foe pie 2 
er «< > 2 yas 3 2 (2.176) 


That implies that the horizontal gradients will be neglected 
with respect to the vertical gradient except in the continuity equa-e 
tion and in the terms involving the horizontal gradient of the mean 
pressure. 


2) The hydrostatic approximation will be used. 


2e9e1 Equations for the Mean Quantities. 


We use the standard coordinate system (x,y,z) in which 2 is 
the vertical coordinate and x and y are the horizontal ones. U and V 
are the mean horizontal velocities and W is the mean vertical velocity. 


The equations for the mean properties become, 
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269-2 Equations for the Turbulent Quantities. 
We define the diffusion operator as, 


D.()=2 ee on (2.183 


KOZ 


and the production components as, 


P, 5 Fe, 24 (2.184) 


The equations for the turbulent quantities can be written 


in the boundary layer as, 
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~wo=— (an? 5 - E69.) (2.197) 
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We now have a closed system of 22 equations in 22 unknowns. 
This set of equations can be reduced to a smaller one by elimination 
of some of the variables. We will also use the turbulent eddy diffue 


Sivities as defined by, 
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- 76 =x, 22 (24203) 
-wq=K, xa (2.204) 


In the case of level 3 the eddy diffusivity for momentum 
has the same value in the x and y directions; this can be proven 
easily by finding the same expression for K, from equation (2.192) 


and from equation (2.194). The division of (2.201) by (2.202) gives, 
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In order to reduce the number of equations, we eliminate some 


variables, use the eddy diffusivities and equations (2.205)-(2.206). 
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2e10 One-Dimensional Level 3 Equations. 


The equations for the unidimensional case are obtained by 
assuming horizontal homogeneity which implies that there is no hori- 
zontal advection. In this case the vertically-integrated continuity 
equation yields, 

W = constant (2.207) 

In fact, because the boundary conditions require that the 
mean vertical wind vanishes at the ground and at the top of the bounda- 
ry layer, the constant W=0. Therefore, the Lagrangian time derivati- 


ve becomes equal to the local time derivative, 


a a om 0 eS 
D >t + U ztvV2o+w eat (2.208) 


In order to simplify things further we will assume that 
the geostrophic wind remains constant throughout the period of time 
considered: this allows us to replace the horizontal pressure gra- 
dient by the geostrophic wind. Therefore, it is not necessary to 
compute “iy mean pressure in a one-dimensional model. We now write 


down all the atmospheric equations that we need, 


2210.1 Eguations for the Mean Quantities. 


oP - 2. (K, SE) =e =v) (2.209) 
— 2. (K, 2H )=£ (U.- 0) (2.210) 
20.2 (x, 38) =o (25241) 
92.9 (x 282) =0 (2.212) 
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2.1062 Eaquations for the Turbulent Quantities, 
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The number of variables in the one-dimensional problem has 
been reduced to 12. We have 11 equations from (2.209) to (2.219) and 
the equation for Q is given by (2.128), which gives us a closed set of 


equations. 
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2ei11 Boundary Conditions. 


In a one-dimensional model there are only two boundaries: the: 
top at which the meteorological variables reach their synoptic values 
and the bottom which is fixed near the surface of the earth. 

2011.1 Upper Boundary Conditions: 220. 

1) All of the turbulent quantities including the coefficients of eddy 
diffusivity vanish at the top of the boundary layer. 

2) All of the mean quantities except pressure are kept fixed at the 

top of the boundary layer. The mean wind reaches there its geostrophic 
value. We need pressure only in the computations of the atmospheric 
infrared radiation which do not require a great precision in the 
computation of pressure so that it is sufficient to determine the 
pressure from the hydrostatic equation. The pressure will be kept fixed 
at the ground which implies that the upper pressure will vary 


with any variation in the mean virtual temperature of the layer. 


Ulzrop) = UZ (2.220) 

Mizees) =V, (2.221) 

Olz op) = synoptic value (2.222) 

Q( zt op) = synoptic value 25223) 

P( 24 op) = value determined by the hydrostatic equation (2.224) 
Turbulent quantities (24 op) = 0 (2.225) 


21122 Lower Boundary Conditions: 2-70. 
The basic assumption is that turbulent diffusion is dominant 


in the surface layer so that the equations for the mean quantities 


simplify to, 
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£ (vw)=0 Vv w = constant = - ue sino (2.227) 
a (we )=0 we = constant = - 9 (2.228) 
ae (wq)=0 Wq = constant =-E (2.229) 


where x is the roughness angle made by the friction velocity u,. with 
respect to the x-axis. H is the heat flux and E is the flux of water 
vapor. Therefore, the formation of a region of constant stress and 
constant heat and water vapor fluxes close to the ground is a direct 
consequence of our assumption that the dominant process is diffusion. 
The next thing to consider is the behavior of K,, the eddy diffusity 
for momentum, in the surface layer. Its dimensions are | em2/ sec | = 
speed x length. The scaling length becomes f= ky (2+ 2,) in the 
lower boundary layer and the appropriate scaling velocity is u Meee 
Therefore, it is appropriate to expect that, 

Kal2-72,) =uyd =k, ty (2 + 29) | (2.230) 

Using the definition of KA? we can write, 
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We integrate the last expression and obtain, 
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V can be derived similarly, 
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The last expression was obtained by relating K, to K, with 
the help of the turbulent Prandtl number which is simply the ratio of 
K,, over Ky. The Prandtl number has to be evaluated from the equations 
in the constant-flux layer. We must be very careful in doing so as the 
following discussion will prove. The equations for Ss and Ky are (2.217) 
and (2.218). In the constant-flux layer we can negleet buoyancy and 
the diffusion terms containing ),, so that the limiting value of the 


ratio of K over Ky is, 


3 3 
= = Saget = (1-3C,) = 0.832 (2.235) 


However, when we use the level 4 model of Mellor and Yamada 
(1974) we get a slightly different result. We start from their equa- 
tions (37c), (38b) and (40c) in which buoyancy and diffusion terms have 
been neglected as well as any time derivative. 

za 2 3) 
8 2 e 8 
Q=- ( woe =) -2 ——. (2.236) 
a 3 i 
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o=- (wee, o%) 32 -24* (2,257) 


2 3 Ny 
on @ SO. sue (2.238) 
We isolate the desired covariances, 
we & - 2 » e” (2.239) 
wes 341 (w2= 0,02) 22 (2,240) 
ee eee (2.241) 


The last two equations give the value of the coefficients 


of diffusivity, 
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K. = 3 f2 we (2.243) 


Their ratio becomes, after using the expression for w 


Pine cae 4 ae 1 
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3 iv 


The difference between the two expressions is caused by 
the rearrangement or terms which occurred when we separated the iso~ 
tropic part from the anisotropic one. In level 3, the dissipation 


term in A, is not explicit but is contained in the terms in e, 


Therefore, the most straightforward way of obtaining Pre is to go 
back to the equations of level 4, 
If we integrate (2.234) we obtain, 
H Pry 2Z+Z, 
O-Oe.) "GR (=F) (2.245) 


Due to the fact that the scaling lengths for water vapor 
were chosen equal to those for temperature, the ratio K/K is equal 
to the Prandtl number. We integrate the equation for the gradient of 


Q obtained from the definition of K, and get, 


E Py 2+z 
=| se 2 
Q¥-1G(2" ee K In ( Zo ) (20246) 


After having multiplied (2.185) to (2.183) by / , we note 
that the diffusion terms in ) and the buoyant terms contain J 
which goes to zero as 2 goes to zero. In the constant-flux layer these 
terms are neglected as well as the small time-derivative terms. The 


final lower boundary values ares 
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e*(z,) = B, w, (2.247) 
ao a H® Be 
2 B | 
a2(2,) = oe (2.249) 
1 
a B 
Rah rep (2.250) 


ue By” * 


u2(2,) = e*(z,) ( — > 1 (4 cos2@xe 2 sin’ o ) (2.251) 
— A 
v'(2,) 2 e*(2,) ( nt 5 (4 sin? x - 2 cos’ x yi G2.252) 
om A 
Wag) = (29) (3 > 2 Be) (24253) 
oe Dees 
uv (z,) = 6 sine cos% e (2) By (2.254) 
as A 
u@(z,) =i "24 ( Pat +1) cos am (2.255) 
By 
pea ha A 
ve(z,) = =A! ( Pie +1) sine (2.256) 
1 
dae ork 
Taz.) = Sey B (P,, + 1) cos (2.257) 
1 
anes A 
watz) = — (P,, + 1) sine (2.258) 
i 
K (2 29) # ke (z5*z) Ux (2.259) 
u 
K,(z->2,) = Kk. (z+2z) Pre (2.260) 
K,(22%9) = ky (24a) (2.261) 
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CHAPTER III 
FINITE-DIFFERENCE EQUATIONS FOR THE ATMOSPHERIC VARIABLES 


3.1 General Remarks. 

There are two main methods for solving numerically the atmos- 
pheric equations: the spectral method and the finite-difference method. 
Each method has advantages and disadvantages. These are discussed 
very briefly below, 

The spectral method involves the decomposition of the varia- 
bles into a finite number of spectral components. This can be accome 


plished, for example, by a Fourier time series, 


- keko yet 
u(x,t) HS v5(kyt) exp(i k*x) (3.1) 
k=1 

where u(x,t) is the variable in the space-time domain and, 

v4 (kyt) is the spectral component in the wavenumber-time 

domain, 

ky is the highest wavenumber that is retained. It is a measu- 

re of the resolution of the spectral representation. 

The atmospheric equations are transformed into wavenumber 
Space which gives us an equation for each of the spectral components, 
The first problem arises when we wish to have a vertical resolution of 
less than one meter in the lower boundary layer. We would need to com- 
pute over a thousand spectral components which is an enormous task when 
the atmospheric equations being used are complex. One way of getting 
around that problem would be to divide the boundary layer into a 10- 


to 50-emeters thick surface layer and into an upper layer extending up 
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59 
to the top of the boundary layer. The spectral equations could then be 
solved separately in these two regions, which would increase the verti 
cal resolution to the desired value in the surface layer without having 
to compute more than 10 to 100 spectral components. The linkage between 
these two layers could present some problems related to the fact that 
we are not solving the spectral equations for the same wavenumbers in 
both layers. 

The finite-difference method has been used for a long time 
and numerous schemes have been devised ranging from simple ones to 
sophisticated ones. Basically the method requires replacement of the 
derivatives by finite differences. We will use centred finite diffe- 
rences., For example, the vertical derivative at the grid point k of 
a function A is approximated by the difference between the value of 
the function A at one grid point above k and its value at one grid point 


below, divided by the distance (2AZ) separating the two grid points, 


A = NR _ AUkt1 “Alki 
34\, e AAI 2k TORE (3.2) 


All of the finite-difference formulae are referred to a grid. 
We can choose the grid as we please. We know that some of the atmos- 
pheric variables, for example, the mean wind and the mean temperature, 
vary rapidly in the surface layer and more slowly above that. 
Therefore, we need much more resolution near the ground than at great 
heights. The ideal solution would require logarithmically-spaced 
grid points in the surface layer and linearly-spaced grid points in 
the upper boundary layer. This type of grid would have the advan- 
tages of giving good resolution near the ground and of involving 


relatively few grid points. This can be achieved in a continuous 
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60 
manner by the log-linear transformation of coordinates with properly 
chosen parameters. There are some drawbacks related to the usage of 
finite-difference methods of which non-linear instability, aliasing 
and errors due to space truncation are major ones. 

The finite-difference method was preferred for its simpli- 
city, its potential for increased resolution near the ground, and its 


smaller computing-time requirement. 


362 Log-Linear Atmospheric Grid. 
The best way of designing a log-linear grid is to transform 


the linear coordinates (z) into log-linear coordinates (Y). This has 
the advantage of simplifying the finite-difference expressions which 
would be complex and position-dependent in a nonelinear grid. In the 
transformed 4 -coordjnates the spacing is linear. An appropriate trans- 


formation for meteorological processes in the surface layer is, 


ries 
Y =A (gn ( P+) (3.3) 


te) 


where, AC = an arbitrary constant that depends on the number of 
grid points and the height of the uppermost grid point, 
He. = a constant chosen as closely as possible to the average 
value of pos the maximum size of the eddies, 
Zo = a constant chosen equal to the roughness height Zo 
whenever possible, 
Ky = Von Karman's constant. 


We will now demonstrate that we have the best grid possible 


in the surface layer it}: =f, and Zo = 25+ If we take the derivative 
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of Y with respect to z, we obtain: 


alee i toa tating ab i 
aa Ae ( ko(atet) * oF ) (3.4) 


The mixing length can be expressed as, 


1 
y) = as LER (355) 
k,(z+zq) fin 


Therefore, in the ideal case we haves 


ees (346) 
and expressions like, 
og aoa are (3.7) 


are encountered frequently and evaluated exactly by the finite-diffe- 
rence method in Y -coordinates. Additional proof comes from conside- 


ration of the wind profile in the surface layer: 


eu As ae oe (3.8) 


ko (ztz,) .71 


pL ck acidity Ae yh ) (3.9) 


a4 AC Le 
As long as (2+z,) << dees U will vary almost linearly in 
Y -coordinates from one grid point to the other, Therefore, the finite- 
difference approximation to (3.9) will be very accurate. Similar 
arguments are valid also for V, () and Q. 


The largest errors would be caused by 2bFzos in which case 


(3.9) becomes, 
du _ x cosm (2 + 25) Ko(z+z,) 
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highly non-linear for small values of z. 
The influence of a difference between J and \: is considerae 
bly smaller in the surface layer than the difference between Z, and 25° 


In this case we have, 


( yu y' _ Ux COSK (ois Kg(ztz,) 


oY a cio Ts ) (3.11) 


(s tps MS (2+2z,) 
pe 
Shes fue e Ce -})) (3.12) 


For example, given Za en, = 27 m, and k, = 0.35, then 
a maximum value for ), is about 8000em, in which case (3.12) becomes 
1 + 0.00008588 (z+z,). The error is less than 1% for 2<¢1m, and it 
is less than 10 % for 210m. Therefore, the influence of a differen- 
ce between Le and |) | is more important above 10 m but in this region 
(3.8) begins to break down. Therefore, we expect that the log-linear 
coordinates will give excellent finite-difference approximations to 
the atmospheric equations in the surface layer without requiring too 


many grid points. 
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63 
303 Atmospheric Equations in LogeLinear Coordinates. 
We will transform all of our equations into the log-linear 
coordinates. 


~e) + 2.) 


Zo £! 


In the transformed coordinates Y= 0 at z = 0, which corres- 


Y= ac ( Lin ( 
(s) 


ponds to the roughness height Zoe AC is a constant which is defined 
by the specification of the number of grid points and the height of the 
top of the boundary layer. If we call the last grid point Y, Lae and 


the height of the top of the boundary layer Ztop» we get, 


2 i Zmaxtzo Geax -1 
AC bare ( ko ln sre + am (3.13) 


The next step is to determine the intermediate values of 2 
by an iterative method based on the behavior of the log-linear trans- 
formation at small and large values of - Close to the ground the 
coordinates are predominantly logarithmic. Therefore, as a first 
approximation we neglect the linear term, 

(a = 2, (oxp(ky4 /AC) =1) (3.14) 
The next guesses involve the complete formula, 
2 = 2, { exp [ k,(4Y/AC - 2/ 4°) -1}} (3.15) 
The iteration is completed whenever the difference between 
two successive guesses is smaller than a predetermined value. Near 
the top of the boundary layer the coordinates become almost linear. 
As a first approximation we assume that the grid points are linearly 


spaced, 


= 016 
2 2 op /Y max (3-16) 
The successive guesses are given by the complete formula, 


2 = 25 { Yiac ~ tel an ( 5] (3.17) 
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The iteration is continued until the difference between 
two successive iterations becomes smaller than a predefined value. 
The complete numerical scheme is more complex if we want a rapid solu» 
tion at all grid points and especially at those intermediate grid points 
where the influence of the logarithmic and linear terms is about equal. 
In Appendix A we describe the complete numerical treatment of the 
derivation of the log-linear grid. 

The first derivative of 4 with respect to z has been obtai- 


ned in (3.4) and the second derivative of Y with respect to z is, 


Ne - AC 
a2 = ~ Klatz )% (3.18) 


The second derivative does not involve | ¢ which means that 
the second derivative is evaluated as exactly as possible as long as 
Zo > 2 The two derivatives are computed by inserting the value of 
zZ corresponding to each grid point. 

Now we will transform the equations from the z-coordinates to 
- the Y -coordinates. In order to generalize the transformation we will 
use the variables A and B as representations of the atmospheric varia- 
bles. Any variable which is not differentiated will keep its identity 
unchanged at the grid point 7, eorresponding to the height 249 

A(z,) = A(Y4) (3.19) 

The first derivative in z-coordinates is transformed into 

y-coordinates with the help of the chain rule of differentiation, 


SALSA DY (3.20) 
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We must be careful when we transform the second derivative 
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transformation is, 


= 2 (24) =2. (242%) (3.21) 


The chain rule of differentiation is valid only for the 
first derivative and has been applied to the inside term of (3.21). 


Now we differentiate both terms with respect to z, 


Seeha) Sey Magoty alike 
re oeay et Ty 2 (3.22) 


After a few manipulations the final result is, 
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Now we can estimate the diffusion terms, 


a SE eA Beg ee 
z fA ¥z’ Uz YZ ape (362t) 


If we use (3.20) and (3.23) we obtain, 


2 (a 38 = 38 (34 (2%)? +a XH) Es 3B a (24)7(3.25) 
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We can now proceed with the transformation of coordinates, 


3e3e1 Equations for the Mean Quantities in Log-Linear Coordinates: 
First Form. 
We will transform the set of equations for the mean quanti- 


ties +(2.5209).6 102.212) in Y -coordinates. 
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33.2 Equations for the Mean Quantities in Log-Llinear Coordinates: 
Second Form, 

Equations (3.26) to (3.28) have been used successfully by 
Estoque (1973) in his numerical model of the boundary layer. However, 
he has used an arbitrary formulation for K, which is less critically 
dependent on the sign of the temperature gradient than the formulation 
obtained from second-order modelling. Some of the problems encountered 
with the utilization of (3.28) are: formation of many artificial invere 
Sions and the blocking of the daytime inversion height at a grid Bot 
near the ground. These problems can be overcome only by very arbi- 


trary assumptions. These problems are discussed in detail in Appen- 


dix C. 
There is another way of writing the diffusion term in y= 
coordinates, 
2. (x, 22 [2 98 d4)] 34 
Td an sy ra oY (Ky t oy 9 =) 2 


This formulation uses the vertical variation of the heat 
flux which insures a consistent vertical profile of not only the 
temperature but also of the heat flux. We can use the same type 
of formulation for U, V and Q and obtain a new set of equations for 


the mean quantities. 


ib fr ( Ky 7 )] ak =£(V-V,) (326!) 
ay - (35 (3p 24 )} ai = £(Ug-U) (3.278) | 
2B fig n OH 24 = 0 (3.28") 
29 Le oe eee (3429") 
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323-3 Equations for the Turbulent Quantities in Log-Linear Coordinates, 
We will now transform (2.213) to (2.219) inito\g7 -coordinates. 
We define the diffusion Bear ) as, 


2 Duayeo ay ee 
dc) (of Boh ata L)*+ 9) Sat) (3.30) 


The equations for the turbulent quantities are: 


‘dod TOs LA wear il 34,2 Kk @.., 3 
vt = 5 Die”) = 2[ Ky jav vy oe rare? = ¢ area ert Ky af as 
= Oe 
- or] (3.31) 


2 ~~ 2 AS PD 
— = 0.23 D(02) = 2[ (Sp) (24). aa (3.32) 
> 42 2 q2 
22° ~ 0,23 D (a) = 2[K, (38 <p) “(au)” By] (3433) 


aes - 0.230) (46) = (KK RE S ; De Oe Ly? ae (3.34) 


A aK 1-3C, JePtix0. 23a, 4 D (02)=21A,2(m 20) +0, 61K, ee) 
Si ites aaa) NET Wa LEE 


o@+6a2 |? al" Gie +9a% \2 (qe 3 7 40.6128 3p 2 (3.35) 
K,=[a, 4 (o7+4x0.23A, 2 D (62)-6A, 1 an (24)?.12x0.614, 0 ex, 22 24 


94 32 
~305(0°/T, +0,6195)/(38 24) = (o2+1207 J2g 2 D 40 5) (3436) 


K=[A, 4 (e344x0.23A, 2 D (02)-6:, | Kn $e : ap oy aed er eo 


sox /1,1046:0°)/(38 So) = (oP ar2n] 170.6105 $4) (3437) 
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3.4 Finite-Difference Equations, 
A Taylor's series expansion around the grid point k, at the 


time t Bs is used to derive the finite-difference equations:s 


x blo maaoiad (A'S iat i ¢ : 
f(x,tm Ax) = f(x.) + mAx a + + nf 35n teooe (3,38) 
Xo Xo 
where, X, can represent te or Ko» 


A x is the increment At orAk, 

m is the number of increments (positive or negative), 

n is any positive integer. 

There are many types of finite-difference approximations to 
the derivatives. For example the forward-eEuler finite-difference 


approximation to the first time derivative uses (3.38) with m = 1, 


L(t ts) ate 
= ee Sone ate (3.39) 
° 


The truncation error is estimated from the first neglected 


term in the Taylor's series, 


ce] 


The centered finite difference is another important type 
of finite difference. In this case wo use (3.38) for m=1 and m= -1 


and subtract the two series from one another to get, 


af] _ f(tytdt) = F(to- At) 


ry’ (3.441) 
P) Ete 2At 
In this case the truncation error becomes, 
2 
ot, = (AE? Yt (3.42) 


rune Bs t2 t, 
which is expected to be smaller than (3.40) because it involves a 


higher derivative term, 
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34-1 Vertical Derivatives, 

We will use centered finiteedifference approximations for 
all of the vertical derivatives, mainly because we expect that the 
truncation errors associated with centered finite-difference schemes 
are smaller than those associated with nonecentered schemes. If we 
make use of the fact that Ay= Qk = 1, the simplest centered finitee 


difference approximation to the first ¢/ -derivative becomes, 


| f(k,+1) = £(k,<1) 
sy Lee ee (3.43) 
° 
whose truncation error is, 
i ieee i 4 
Tune * 7 6 34/3 x. (3044) 


Therefore, (3.43) will be in error whenever the third 
derivative is important. We have already shown that the mean atmos- 
pheric variables vary almost linearly in Y -coordinates in the lower 
boundary layer. Therefore, the third derivative is small in this re- 
gion and (3.43) is an excellent approximation to the first derivative. 
However, at some height above z, we expect large non-linear variations 
in the atmospheric variables. In this case the third derivative becoe 
mes relatively large and causes serious truncation errors. One way of 
reducing this effect is to use more grid points which is very costly 
in terms of computing time. An alternative approach is to use a higher- 
order approximation to the derivative in the region where the non- 
linear variations are expected to be important. This is done by adding 
the finite-difference expression for the truncation error (3.44) to 
the finite-difference expression for the first derivative (3.43). 


When this is done we obtain, 
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Br[ oo = tlegta) + 8 [e(egtt) ~ £(koe1)) + £(ko-2) aes 
94 1k, a 


and the truncation error becomes, 
cise ome (3.46) 


The atmospheric variables vary relatively smoothly in the 
boundary layer and the fifth derivative can be expected to be fairly 
small. This can be checked by evaluating the fifth berivelies by a 
finite-difference expression. 

In order to find the second derivative we form Taylor's 
series for m=1 and for m= <1, add the two series together and 


obtain, 
Qa = f(k +1) -2 f(k.) + f(k-1) (3.47) 
34 2 ° ° oe : 


The truncation error is, 


4 
k = 1 f 
setae) camry Th (3.48) 


The truncation error of (3.47) seems to be fairly small be- 
cause it depends on the fourth derivative. However, in the case of a 
strong nocturnal inversion (3.47) may have a relatively large trunca- 
tion error at a few grid points. If this is needed we can use the 


following improved approximation to the second derivative, 


32 f) «of (kg +2) 416£(k ott )=30f( kg) +16f(kg-1)-f(Kp-2) Guia) 
94 Zl Ko 12 ; 


whose truncation error is, 


6 
k Peet Oat 
Tune erg 346 (3.50) 
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3.4.2 Time Derivative. 


Suppose that we must solve, 

Sz = A(f,t) (3.51) 

We estimate this expression at the time t,+$ At. The right- 
hand side of (3.51) can be estimated at the time t,+3 At as the ave- 
rage betwoen the value of A at ty, and its value at t,+At, 

A(f,t +2 At) = S[A(f,t,) + A(t t th t)} (3.52) 


and the leftehand side becomes, 


Qf _ f(tot At)=f(tg) 


Equation (3.53) is similar to (3.39) but they are taken to 


be valid at different times. The truncation errer becomes, 


2 3 
th Oboee SUA Rt oe , 
Trune ahs G3 tte at (3.54) 


This mathod is implicit because it involves the function 
A at the time step t.tAt which is not know a priori but must be 


obtained by some iteration technique, 


304.3 Stability of the Finite-Difference Equations. 
Many of the equations in the model are of the forn, 


of =a) EE £+a(2) 2£4+4(3) £ (3055) 
Re ay% OY 


where A(1), A(2) and A(3) are functions of f, t, Y and of other 
atmospheric variables. The equations are, therefore, non linear 
which renders very difficult a complete analysis of stability. We 
will assume that the coefficients are constant in order to simplify 


the analysis of the stability of the finite-difference equations, 
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We rewrite (3.55) using the finite-difference approximations as given 
by (3255) (3.52) (3.47) and (3.43). 
£(ko st ,+At)-£(kost,) 


SSCA Es em A(1) [fUkgtt yt )=2f (Kost y)tf(koe1ety tf (Kot sty 
+At)e2f(kostot A t)tf(koelytot A t) }4A(2)[f(kottsto)=f (komt sty) +f(Kotts 
tri-£(k -1,t + At)] /24A(3) [£(k, st, )+f(KostotA t)] (3.56) 

We suppose that the exact solution of the finite-difference 


equations is one of the Fourier modes, 


f(kotm,t,im At) =D exp\i[(kytm)l-(toin at)o}} (3.57) 


where, D = amplitude, 


f 


1 = wavenumber, 
= frequency. 


Substituting (3.57) in (3.55) and simplifying, 


eat sgh 2 bay) 401) (ctttemtt ny 4a(2)(ottno tt) 44(3)/2 
(a t/2) (3.58) 


We replace the exponentials by the appropriate sine and cosi- 
ne terms and obtain, after sorting out real and imaginary parts, 
cos cht -1=At[cos(oAt) A(1)(eos 1~1)#A(3)/2\=a(2) (sin vat sinl] 
(3.59) 
sincb t=ht|sin(s At) A(1)(eos 1 -1)+A(3)/2]+A(2) coscdt sin l 
(3.60) 
We square (3.60) and obtain, after a few manipulations, 
1 


sin* cAt = {- At/A(1)(cos 1 @1)+A 22\" 
1+ CoACS eSsineL 


Therefore, the following is true for any real l, 


sin®gAt< 1 (3.62) 
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This implies that ois always real and that the computations 


are unconditionally stabdle. 


34.4 Comments on the Finite-Difference Equations. 
We need initial conditions in order to use the time-dependent 


model. If we were devising an operational model we would need initial 
conditions that correspond closely to physical reality and we would 
have to rely on measured atmospheric variables. In our case we do not 
have any experimental data for initialization of all variables. We 

are interested in a numerical simulation which would represent a possie 
ble urban heat island but not an actual one. Therefore, it seems justi- 
fiable to use the steady-state equations in order to derive the ini-~ 
tial conditions. However, if we use the steady-state equations for 

all atmospheric variables, we are led to an undesirable result. The 


steady-state equation for Gis, 


Slay ago hw 

aia UNE ae (3.63) 
which implies that, 

Ke ) = constant (3.64) 


Therefore, the assumption of steady state implies that the 
constant-flux layer extends up to the top of the boundary layer, which 
is physically unrealistic because we do not expect the formation of a 
constant-flux layer deeper than about 100 m. Moreover, (3.64) implies 


that, 


oD . constant (3.65) 
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This would cause the temperature gradient to increase near 
the top of the boundary layer where K, decreases to a very small value, 
One way of getting around that problem is to assign a thickness to the 
constant-flux layer and to keep the temperature profile fixed above 
that layer. This has the disadvantage of being rather arbitrary, An 
alternative is to assume an isentropic atmosphere, because we expect 
that the atmosphere will be nearly isentropic at the time of maximun 
temperature. 

The equation for Q is similar to the equation for @) ° there: 
fore, we suppose a constant mixing ratio initially. Evaporation de- 
pends on the temperature of the air and it is very likely that the 
maximum value of Q occurs at the time of maximum temperature ifthe water 
supply is not limited. Therefore, it seems reasonable to expect that 
the latent heat flux and the sensible heat flux will behave similarly, 

The above discussion is valid certainly for the lower 
boundary layer where the changes in the surface properties are rapidly 
felt at all levels. However, above 100 m there is a time lag involved 
so that at the time when the heat flux vanishes at the surface (near 
the time of maximum temperature), the heat flux might still be impor- 
tant above the surface layer. This points out that the initial condi- 
tions can be fairly poor above the surface layer and that the model 
should be run for a few days of real time before we are certain that 
the model is independent of the initial conditions, 

The equations for the steady state are considerably simpli- 
fied by our assumptions and only U, V, Kw e and { need to be computed. 
The transformation of all the equations into finite-difference form is 


a fairly lengthy exercise which adds no physical insight to the pro- 
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(é) 
biem. For this reason the transformations are presented in detail in 
Appendix B for the steadyestate equations, and in Appendix C for the 


time-dependent equations. 


305 Steady-State Model. 


We will use the Gauss-Seidel iteration scheme. Let Feomp? 


Fast? and Fyows be the value of F calculated from the equations, the 


last value of F in the iteration scheme, and the new value of F in the 


iteration scheme, respectively. We call the difference, 


DINE, (=ak hat cane (3.66) 
and let the new value be, 
|e cates Foast + OMEGA DIFF (3.67) 


where OMEGA is the acceleration factor which can vary from 0 to 2, 
OMEGA = 0 we keep the last values F... = Fast 


OMEGA = 1 we'keep the calculated value FL. = Foast 
OMEGA >1 acceleration 
OMEGA <1 deceleration. 

The iteration needs some initial values in order to start 
the computations. It is preferable to use the best possible estimate 
of initial conditions. For example, we can suppose a linear variation 
of U and V from 2, to 1000 m, where they would reach their geostrophic 
valuese Above 1000 m U and V can be assumed initially to have their 
respective geostrophic values. We can assume a constant initial value 
of about 50000 om”/sec for Ka in the lowest 1000 m and a linear de- 


crease from that value to its top boundary condition above 1000 m. 


We can assume a value of around 3000 cm ror }, and computers from its 
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76 
equation. We initialize 6 by assuming a linear decrease from a maxi-e 
mum of the order of 80 em/sec at the roughness height to its minimun 
top boundary value. 

In the lower boundary layer and up to about 1000 m an OMHGA 
of 1 to 1.2 is satisfactory. When we use (3.26) and (3.27), we find 
that OMHGA must be reduced in the upper boundary layer if we want cone 
vergence of the iterations. K, is determined partly by the vertical 
gradient of the wind and the wind is determined partly by the verti- 
cal gradient of K,. Near the top of the boundary layer the wind vie 
ries very slowly with height so that great accuracy is needed in the 
wind in order to estimate properly K,. For the first few iterations 
we can expect large errors in the vertical derivative of the wind 
and these errors can cause K, to be estimated incorrectly by up to 
one order of magnitude. This in turn increases the error in the wind, 
and the cycle is repeated over and over and causes rapid divergence. 
Therefore, in order to get convergence we have to use a small OMGA 
near the top of the boundary layer. However, when we use (3.26) and 
(3.27) these difficulties do not arise because these equations com- 
pute the stress directly and, therefore, will give reasonable values 


for the vertical gradient of the wind even if K, is badly initialized. 


3.6 Time-Dependent Model. 


We use an implicit scheme which we have shown to be uncone 
ditionably stable for linear equations. A time step of 10 minutes was 
used succesfully and is a fair compromise between resolution and 


computing time. It would not have been possible to use explicit me- 
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77 
thods which must obey criteria limiting the time step because the ver- 
tical spacing in the lower boundary layer is very small and would have 
required extremely small time steps. The Dufort-Frankel (1953) scheme 
was tried unsuccessfully because it showed instability by producing 
waves of large amplitude and of period 2At. These waves could not 
be eliminated even after using a time filter which removed waves of 
period 2 At exactly. This instability was rather small when a time 
step of 1 second was used, but such a time step required too much 
computing time. 

The variables at the time t, +At were evaluated by an ite- 
rative scheme which used the accelerated Gauss-Seidel method, 
The equations were relatively well behaved in the surface layer where 
an OMEGA slightly greater than 1 could be used generally. However, 
during the nightime inversion K, decreases rapidly in the region bet. 
ween 5 mand 100 m. This causes some problems in the iteration scheme 
similar to the ones encountered in the solution of the steady=state eu 
quations near the top of the boundary layer. Therefore, an OMEGA of 
about 0.7 has given good results whereas an OMEGA greater than 0.8 has 
often caused divergence of the solution. Near the the top of the boun- 
dary layer the variables converged more rapidly in the time-dependent 
model than in the steady-state case because of the domination of the 


time-derivative term which tended to stabilize the iteration procedure, 


The best value for OMEGA has to be found by trial and error. 
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307 Initialization at Each Time Step. 

In order to reduce computing time it is desirable to use 
the best possible initial values for all the veriables. This can be 
done by extrapolating from the values of the variables at the past 
time steps. If we stored only two time steps we must use linear 
extrapolation which could be inadequate in some cases of rapid chan- 
ges as happens near sunrise, for example. A quadratic extrapolation 
would certainly give much better results et the expense of an in- 
erease of 50% in the storage requirements. We now derive 2 forma 
for the extrapolation. The variable A can be expressed in terms of 
a second degree polynomial like, 

A=at®+btte (3.68) 


We compute A at t = ©, 1 and 2, 


a(O) =e (3.69) 
A(i)=at+bte (3.70) 
A(2)=4at2b+e (3.71) 


After elimination we find that, 


- = A(2) - 2 AG) + af) (3.71) 


2 
& AC i) = 3 A(O) = Al2 (3.72) 


The extrapolated value at t = 3 becomes, 
A(3) =9atZbte (3.73) 
We use the values of a, b and ec in (3.73), 
a(3) = a(o) + 3 [a(2)-a(1)] (3.74) 
Formla (3.74) is used to get the extrapolated value of 


all of the variables. 
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CHAPTER IV 
SOIL LAYER 


4,1 Introduction. 
The soil layer is modelled relatively simply. The basic 
equation is, 


DT; 9 o 1s 
Dt ae 3 oe (K, 5) my (4.1) 


where, T; = soil temperature, 
Kg = soil thermal diffusivity, 
Zz. = vertical coordinate into the soil. 


Equation (4.1) simply implies that the changes in soil tem 
perature are caused by the thermal diffusion between adjacent layers 
of soil. The basic assumptions are, 

1) The soil composition is vertically homogeneous so that K, is 
constant for the soil layer. However this assumption can be 
removed easily in the case of a more complex soil layer, as 
Ks can be varied in any manner without basic alteration of 
the model. For example, we can superpose layers of soil 
with different iWoemnleds eareiea tyesy or we can vary K, with 
the help of any function representative of the vertical varia- 
tion of K, in a given soil layer. 

2) ‘The horizontal thermal advection is negligible. This implies 
that the mean thermal gradient is vertical and that the soil 
layer is fairly homogeneous in the horizontal at all levels. 


3) The lower boundary is fixed at 1 m below the surface of the 
79 
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80 
earth where the temperature is assumed to be constant at all 


times over the period of one day. 


4,2 Logarithmic Grid for the Soil Layer. 
We expect that we will find the steepest temperature gradient 


in the top few centimeters of the soil layer and the weakest gradient 
near the lower boundary. It is therefore appropriate to devise a grid 
which gives more resolution close to the earth's surface than at im 
below it. A logarithmic grid has been devised for the soil layer ea 
makes use of the following transformation of coordinates, 
Se =A‘ In (1+z,) (4,2) 
where, €= logarithmic coordinate, 
A‘= a constant determined from the choice of the number of 
grid points and the depth of the soil layer. 
We have chosen a grid with 10 points. The transformation of 
coordinates is completely definsd if we require that, 
S = 0 at 2, = 0 cm 
c = 9 at z, =100 cm 
In this case we have, 
A® = 1.95 


= 0.000000 em at & = 0 


(3 
t 


= 0.669945 em at © = 1 


iS) 
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= 1.788719 cm at one 2 


iS) 
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= 3.657010 cm at & = 3 
= 6.776953 em at € = 4 
11.98709 em at € = 5 
= 20.68774 om at € = 6 
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wS 
t 


= 35221734 em at Es 7 
= 59.48100 em at § = 8 


ss 
} 


100.0000 cm at& =9 


NS 
uy] 


The first and second derivatives of € with respect to z, can 


be evaluated exactly, 


oe 


Zs 


(4.3) 


+b 
a = 


) 

De) = Coxe 

0 2 (1+z,)* ce 
8 


Equation (4.1) is transformed into eercominetes and becomes, 


T. eye (ptatsrae hts, 28,’ ] Tees 


t BS Me aes 


4.3 Equation for the Steady State, 


The equation for the steady state is obtained by setting 


the time derivative to zero, in which case (4.1) reduces to, 


27 
Repos oes o (446) 
d 22 
8 
This implies that, 
at 
© = eonstant = C (4.7) 
0 2, 


Under steadyestate conditions we will obtain a linear 


variation in the soil temperature. 


4,4 Time-Dependent Equation. 


We will use an implicit finite-difference approximation to 
the time-dependent equation so that the computations will be absolute. 
ly stable. The implicit scheme is the same as the one for the atmos~ 


pheriec equations in which the time derivative is evaluated at the time 
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82 
t,t At. We transform (4.5) into a finite-difference equation and 
maka use of the fact that AS = Ak = 1 to get, 
Ty(Kostot At) (1+ At K,DDZ2(ky)]=Ty(Kosty)+ St Kf [Py(Kottsty) Tg (kytt, 
tot At )-Te(komtsty)“Tg(Kon-tstotA t)] DD2Z(ky)+[Te(kotl st, +g (kotlstyt 


At)=21 4 (kosto)=2%e(Kostot At)4g(Kon1sto)Mfg(KgrAytot At)] DDZ2 (Kk, )} 


(4.8) 
2 
where, DDZ2(k,) = ( 28 aE 
fe) 
2 
DD22(k,) = 2 ye Me. 
Zs fe) 


4.5 Initial Soil Temperature Profile. 
The urban heat island model is started at the time of maximum 


temperature because the atmospheric heat flux is approximately zero near 
the ground at that time, However, it is probable that the soil heat 
flux reaches its maximum value near the time of maximum surface tempe- 
rature so that a special initialization is needed for the soil tempe- 
rature profile. This can be done by integrating the soil equation for 
a few days of real time, letting the surface temperature vary according 
to a diurnal cycle as given by a cyclic function like a sine. At the 
end of a few days of real time the soil temperature profile should have 
adjusted to a value fairly close to the one it should have over a 
diurnal cycle. The soil temperature profile computed at the time of 
maximum temperature will be used as the initial soil temperature pro- 
file in the urban heat island model. However, the exact amplitude of 
the diurnal temperature wave is not generally known a priori so that 


the derived temperature profile must be put into a form which would 
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83 
be useable for any temperature difference between the surface and the 
lower boundary. This is done by computing ZS(k) which gives a eine 
to DIFF, the temperature difference between the surface and the lower 
boundary, at each grid point. 

DIFF = 7,(0) = 1,(9) (4.9) 
And 2S(k) is defined by, 
28(k) =|T,.(k) = Tg(9)] /DIFF (4.10) 
Therefore, in the complete program for the steady state, we 
will not compute the temperature profile but assume it to be of ine 
form, 
Ty(k) = 7,(9) + 28(k)x [T,(0) = T.(9)] (4.11) 
Equation (4.11) should give better initial conditions than 


(4.7) which assumes a linear variation of temperature in the soil. 
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CHAPTER V 
SURFACE TEMPERATURES 


The temperature at the earth's surface is determined by 
a heat-balance equation which includes long-wave and short-wave radia- 
tion, heat flux from the ground, latent and sensible heat fluxes from 
the atmosphere and artificial heat input generated by a city. 
R,=- (PC, H +P, cg S +P LE) +A (5.1) 
where, R, = net radiation at the earth's surface, 
P= atmospheric air density, 


Ps 


C. = heat capacity at constant pressure of the atmosphere, 


density of the soil, 


C, = heat capacity of the soil, 


latent heat of vaporization, 


sensible heat flux from the atmosphere, 
= latent heat flux from the atmosphere, 


sensible heat flux from the soil, 


= SCs ho eee 
t 


artificial heat. 


51 Net Radiation. 


The radiation balance at the earth's surface is, 
Roe Lot Ro Rt (5-2) 
where, I, = net solar radiation, 


R |= infrared sky radiation, 


R f= terrestrial outgoing infrared radiation 
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5e1o1 Outgoing Terrestrial Infrared Radiation. 


We will assume that the earth's surface radiates like a gray 
body, 
Rf=€, or (5.3) 
where, € {= emissivity of the earth's surface 
0 = Stefan-Boltzmann constant 


T = temperature of the earth's surface 


501.2 Sky Infrared Radiation. 
Most of the sky infrared radiation is emitted by water va- 


por and carbon dioxide. The exact radiation transfer equations become 
extremely complex when we divide the spectrum into many frequency 
bands in order to model the selective absorptivity and emissivity of 
HA0 and C05» In order to keep the model simple it is necessary to 
average the infrared radiation over all of the spectrum. Thea climato- 
logical models which use only the observed temperature and humidity 
at the observing weather station are too simplistic. A compromise 
between accuracy and reasonable computing time is the Brooks (1950) 
model as described by Atwater (1966). Brooks has presented values 
of the emissivity € and of — » the slope of the emissivity curve 
for water vapor with respect to the path length w. Haltiner and 
Martin (1957) have modelled the infrared radiation due to the atmos- 
pheric carbon dioxide in terms of the infrared emission of the earth's 
surface, 

R} (COz) = 0.1864 0 i! (5.4) 


We will derive the radiative transfer equation for water 
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86 
vapor. Our starting point is Schwarzschild's equation for monochromae 
tic radiation, 

dIy=- Pk, (I, - By (t)] as (5.5) 
where, I.) = intensity of radiation at frequency Vv, 
ky= absorption coefficient of water vapor in a thin layer ds, 
By = blackbody radiance as given by Planck's law, 
Now we relate ds to the absorber path length du, 
ds = sec W du (526) 
where WwW = zenith angle. 
We inelude the integral over the zenith angle in du and 
obtain for the downward radiation, 
al = ky [Iy~ 8)(7)) au (57) 
We integrate (5.7) over frequency and path length, 
o 
PY (u ped = \ oh B Poy) du (5.8) 
where wy = total path length. 


The Brooks model transforms (5.8) into a finite-difference 


equation, 
F,, (0) = (ort), (3 £) Au (5.9) 


wnere Au is the path length in a given layer as given by 


(as fp eG 
Au Cy 1000 q (5.10) 
where, a = constant between + and i, 
g = gravity, 
p = pressure in millibars, 


= mixing ratio for water vapor. 
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Elliot and Stevens (1961) have derived equations for the 
emissivity from the data of Brooks (1950), and Hales et al. (1963) 
have presented an equation for short path lengths. In Table 2 we sumna- 
rize these equations, The constant a which appears in (5.10) has 
to be determined from radiation experiments. Atwater (1967) points 
out that a value of $ or 1 is used generally, In our model we will 
choose the simplest formulation (a=1) because the actual formulation 
for the infrared radiation is not a crucial factor in understanding 
the urban heat island when the radiative flux divergence is neglected. 

In Appendix D we deseribe in detail how the computation of 
the infrared is done. In that computation we have made the assumpe 
tion that the surface pressure is constant over the period of one 
day. This is justified because the actual value of the pressure is 
not used in the computation of the mean wind but only in the compue 
tation of the infrared flux. However, this is not a serious restrice 
tion because we can allow the surface pressure to vary according to 


the synoptic conditions if these are known. 


501.3 Incoming Solar Radiation. 
The solar radiation received at the earth's surface is, 


I, =R, (ica) cos (Z) (ty) (501) 


where, R, = solar constant, 
a = albedo of the earth's surface, 
te = transmissivity of the cloudless atmosphere, 
Z = zenith solar angle determined by, 
cos Z = sin (P) sin (S) = cos (0) cos (5) cos (2USlt) (5.12) 
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TABLE 2. Equations used for emissivity in the Brooks method, as a 
function of total path length (cm). This table is a repro. 


duction of Atwater's (1967) Table 2, p 827. 


E = 0.1579 log(1 + 4275 w) o<w<107t 
E = 0.0396 In w + 0.389 107 *< w<1079 
E = 0.0565 In w + 0.506 1079 < w <1077 
E= 0.0653 In w + 0.546 107?< w <107! 
€= 0.0778 In w + 0.575 107*k w 


The first entry is from Hales et al. (1963) and the last four 


entries are from Elliot and Stevens (1961). 
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89 
where, ic = latitude, 
§ = solar declination which depends on the time of year and 
varies from 23° (summer) to -23° (winter) 


Sl. = earth's angular velocity. 


5.2 Sensible Heat Flux from the Atmosphere. 


The atmospheric heat flux is, 
lim p) 
iS 2-7 Zoe (Ky ae (515) 


This theory is valid only above the roughness height whereas 
the heat-balance equation is applied at the interface between the ate 
mosphere and the soil. Therefore, we need an equation relating the 
temperature at Ze to the temperature at the ground level. There are 
many possible equations depending mainly on the nature of the rough» 
ness elements. The most common assumption is to assume an isothermal 
layer from the earth's surface up to z. This has the advantage of 
simplicity and is a justifiable assumption in the case of small roughe 
ness elements. However, in the case of large roughness elements like 
a forest or a city such an assumption is certainly not realistic for 
that layer of air. In the case of a thick forest wea can expect that 
the maximum temperature will be reached at about the level of the 
average height of the trees. In that case the effective surface 
where the radiation balance should be used would not be at the ground 
but at the height where the vegetation absorbs most of the incoming 
solar radiation. Therefore, we can expect that the modelling of the 
layer between the roughness height and the earth's surface will depend 


mainly on the physical situation to be modelled. 


tienes to net eee ae 


(ERes2 


ve 


a oohet maordaies aie wie Wl er Se 
vewtesi ‘nemeshihak ih Si pe rhagh Ak, need 

nid jecaike iiet tact ne ene ie oe” deal 
ee ee te styec gh a 
yee 20) Tee wndan vx sw tts gel baa 

oak Oe, OS eR, ae ut ssLoegene ae pha oP 4 
> spedaretian nibbvetend aM cas mene it na 

spied tient te eee aed i en ghee Nye cee 
sig t eimyepde wtengdot egiel/ Ss sna ‘atd se pean dale peers 
<bt srfetioer don pitstedted uh on batsman aH gira * ‘Tain aw Sera s 
fead Siang man ob tie inkid 8 bee od nt site "io aged 4008 
oes to Inve! atv - ss te Secinags ne Lin v= tinea? pata AE 

- geatnae entero). ait, ah, Fat at - sascost. 96 3 Aland abetae © 

ftbehy bia do 2S Pew lea tome: ad digesta onan lao pats than edt wxete | 4 


nj? 


on 


‘a 


ae! 


@ 


- 


eaterced rt, ie vad eth o eelthvagpy. ad Rok eee od oe 
: ‘ 


=e 


att’ te yetiinbAe #6) 2) See ABR, peters — talon 
Pisahiecevt Llp keer» si wea iee geet ch a s-mulgent, Pe et 
| ; 

pe ee ems” i025 co £ ec aey oft we tlclen ; 


& 


on 7 


90 
Consider an extreme case in which the assumption of constant 
heat flux would be valid down to the ground. This could be the case 
for a relatively rigid, tall and sparse grass through which the air 
eould circulate easily and which has a relatively high roughness height. 


In this case we can apply (5.13) in the layer below z.. We average 


0° 


(5.13) over the layer between z, and the ground and obtain, 


H =k 28 (5.14) 


In a finite-difference form this is expressed as, 


O(2,) =@ (2=0) +2, ie (5.15) 


The coefficient of eddy diffusivity for temperature K;y varies 
linearly near the roughness height so that a good average for K, over 


the layer is, 


— Kk 2% 
Ke = $[ Ky (z9) + K(0)] = (Kuen rane (5616) 


where, (K,) = molecular diffusivity = 0.18 em“/see. 


mol 
In general we can neglect the molecular diffusivity with 
respect to the turbulent diffusivity, in which case (5.15) ean be 


rewritten as, 


P 
@(s,) =OQx-0) +t w (5417) 


We can generalize (5.17). H is computed near the roughness 
height and is an interesting variable to keep in the formula as well 
as the ground temperature and therefore the following can be used, 

z,) =O(e=0) + BH (5.18) 


Some special cases of this formula are, 
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B = 0 implies an isothermal layer between Zo and the ground, 
2 Pr, 
B= z implies a constant heat-flux layer between Z, and the ground. 
° 


B<O implies that the temperature gradient in the layer between Ze and 
the ground has a different value from the temperature gradient in the 
layer above Zoe This could be the case, for example, for a dense fo- 


rest, 


53 Sensible Heat Flux from the Soil. 
The sensible heat flux from the soil is simply, 


alain ats ) 


S = 750 (K, CTE (5.19) 


5+ Latent Heat Flux from the Atmosphere. 
The latent heat flux from the atmosphere is, 


_ lin 


B= 3, (ky $2) na (5:20) 


5.5 Artificial Heat Genoration. 

Artificial heat generation is likely to be more important 
during the winter months than during the summer months. In fact, 
during the summertime we would expect that the man-generated heat would 
be insignificant when compared to the high solar energy input. All of 
the numerical simulations will be done for summertime situations in 


which artificial heat is neglected completely. 


5.6 Solution of the Surface Temperature Equation. 


We now have a mathematical expression for each of the terms 


ain (5.1) and the unimown variable is surface temperature. The next 
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92 
step is to write down these expressions in finiteedifference forn. 
The resulting fourthedegree equation is then solved exactly yielding 
directly the surface temperature. The numerical solution is described 


in detail in Appendix E, 


ob om D adalk es 
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CHAPTER VI 


NUMERICAL SIMULATION OF THE URBAN HEAT ISLAND 


6.1 Steady-State Model, 
The steady-state model is used mainly to provide initial 


conditions far the time-dependent model. Therefore, we do not require 
a great accuracy in the numerical results because that does not affect 
greatly the time-dependent model. This implies that we must be Rarer 
in our interpretation of small differences between two steady-state 
simulations. Three situations were simulateds 
a) small roughness height (1iem) and moderate geostrophic wind (10 m/sec). 
b) Small roughness height (1 cm) and strong geostrophic wind (20 m/sec). 
ec) Large roughness height (100 em) and moderate geostophic wind 
(10 m/sec). 

From these three cases we should be able to assess the effects 


of the geostrophic wind speed and of the roughness height. 


6.1.1 Effect of the Geostrophic Wind Speed. 
Fig. 6 represents the wind hodographs for the three cases. 


The hodographs for moderate and strong geostrophic winds are virtually 
on the same normalized curve. The main difference is the height at 
which the relationship between U and V is realized. Table 3 shows more 
clearly the effect of the increased wind speed. We define arbitrarily 
the constant-stress layer as the layer close to the ground in which the 
stress does not vary by more than 20%. The constant-stress layer doubles 
approximately with a geostrophic wind twice as strong to reach a thickness 
93 
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Figure 6, Normalized wind hodographe The various symbols are: 


© for the case V, = 10 m/sec, 2, = i em, 


fe] 


A for the case V, = 20 m/sec, 2, = 1 em and, 


GO for the case Ve = 10 m/sec, Zo = 100 cme 
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U/Vg 


The number inside the parentheses is the height in meters, 


Table 3. Comparison between three steady-state simulations. 


case 
parameter #1 #2 #3 
Bo 1 om 1 cm 100 cm 
V 10 m/sec 20 m/sec 


10 m/sec 


26.149 


17 653° 
1.0504 


170349 


1.0630 
0.264 


1.0477 


max U/V, 


0.210 
56.76 m/sec 
240 m 


max V/V 0.205 


30.19 m/sec 45,26 m/sec 
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of 240 m. The roughness angle & is comparable in both cases. Greater 
accuracy in the numerical results would have been needed in order to 
assess with certainty the effect of the increased geostrophic wind on 
Oo. From Estoque’s (1973) numerical model we would have expected a 
decrease of 1° ina due to the doubling of the geostrophic wind and not 
an increase of 0.2° as obtained in our numerical simulation. The maxi- 
mum normalized values for both U and V are larger by about 1% for the 
ease of strong geostrophic wind. 

In Fig. 7 we compare the behavior of e¢/u% for all the situa- 
tions. We observe that the turbulent energy is nearly constant in 
the constant stress-layer as indicated by the lower straight line, but 
decreases rather “pedis: above that layer for the case of moderate geo~ 
strophic wind and much more slowly for the case of strong geostrophic 
wind. Fig. 8 gives the vertical profile of Ky/k, uy. The lower boun- 
dary value is z,. In the constant-stress layer we have a linear rela- 
tionship between K, and z as indicated by the lower straight line in 
the nearly log-log graph. The vertical axis in Figse 7 and 8 is the 
Y -coordinate expressed in terms of ordinary heights: therefore, it is 
logarithmic in the lowest part and nearly linear in the upper part. 
Both Figs. 7 and 8 indicate that the turbulence will be greatly increae 
sed above the constant-flux layer in the case of strong geostrophic 
wind. This causes a substantial increase in the boundary-layer depth 
which exceeds 1900 m for the strong wind case as compared to 1300 m 


for the moderate wind case. 


6.1.2 Effect of Increased Roughness Height. 


In Fige 6 we observe that the rougher terrain causes an 
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Figure 7. Vertical variation of e%/uz t a) over Zo=1 cme Solid lines 


yee m/sec. Dotted lines V_=20m/sece ») Over 2 =100 cme 
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Figure 8, Vertical variation of K,/k, ux t a) Over z,=1 em. Dotted line 


for V=20 m/sec and solid line for ¥,=10 m/sece b) Over z2,=100 cm. 
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97 
increase in the roughness angle as expected from theoretical conside~ 
rations. The increase of 8.5° for a change in roughness height of 2 
orders of magnitude agrees very well with Estoque's (1973) Fig. 6. 
From Table 3 we can see that the maximum values for U and V are, 
respectively, 1.5% and 30% larger over rougher terrain. The top of the 
constant-stress layer is increased by about 50% by a change in roughness 
height of two orders of magnitude. The friction velocity uyis also 
increased by 50%. This is an interesting result because Table 3 seems 
to indicate a linear relationship between u, and the top of the cons- 
tantestress layer under neutral conditions. However, such a relatione 
ship would need to be investigated over a wider range of roughness 
heights and of geostrophic wind speeds. From the knowledge of u, we 


could estimate the height of the constantestress layer by a formula 


likes 
Top of c.-s. layer = 400 uy, 
From Figs. 7 and 8 we deduce that both e* and Kn are increae 
aA by an increase in z,. Again, the increase in the turbulence is 


greater above the constant-stress layer than below it. This increases 


the depth of the boundary layer by at least 50%. 


6.2 Time-Dependent Model. 


The time-dependent model was run with the second form of 
the equations for the mean quantities. The model ran smoothly with 
OMEGA = 0.7. We have supposed that convergence was reached whenever 
the differences between two successsive iterative values for the dif- 
fusivity coefficients were all smaller than 10 em“/sec. In tha cons~ 


tant-stress layer, this generally insures an accuracy of 0.1 m/see in 
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the mean winds and of 0.001°C in q). Above the constantestress layer 
the accuracy is generally greater because the variables are varying | 
more Slowly with respect to time. The model is started at the time of 
maximum temperature which occurs near 2 PM. The parameters common to 
all simulations are: 

b = 23° (sumer solstice) 

Q = 53° (Edmonton's latitude) 

T= 0.85 

E = 0.82 

albedo = 0.25 
P (z,) = 940 mb (altitude near 2200 feat) 

Q = 0.006 g/g 

Four simulations were made: 
a) Moderate wind, small roughness height and high soil conductivity. 
b) Strong wind, small roughness height and high soil conductivity. 
c) Moderate wind, large roughness height and high soil conductivity. 
d) Moderate wind, small roughness height and low soil conductivity. 

Table 4 lists the differences in the external parameters for 
the four cases under study. From these four simulations we should be 
able to assess the effect of the geostrophic wind speed, the effect of 
the roughness height and the effect of the soil conductivity. Using 
simulation a) as a reference, we will compare U, V, @, K, and wo 
for the three other cases. Then, we will examine the diurnal cycle 
of the other turbulent quantities using model b). Model b) is chosen 
in this case because its diurnal cycle is smoother than the others 


due to the increased turbulence caused by the stronger geostrophic 


wind. 
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Table 4, Parameters for the time-dependent simulations. 


parameter 


29 [om] 
V, {m/sec} 
Ps [e/en?] 
k, [em?/sec] 0.021 0.021 


ce fergs/(g- °c)} | 1.25x10? 1.25x10? 


k. ¢ 
Ps Ks %s 4, 53x10° 4. 53x10" 


ergs/(cmesec- °¢)) 
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6.2.1 Geostrophic Wind Speed. 

Fig. 9 shows the diurnal variation of the surface temperatu- 
ree The increased geostrophic wind reduces the amplitude of the diurnal 
eycles the minimum temperature is increased by 2.92 °C and the second- 
day maximum temperature is decreased by 1-71 °C. There is no detectable 
time lag between the two temperature cycles and they both reach their 
minimum near 4,30 AM and their maximum near 2.30 PMe The surface 
cooling rate illustrated in Fig. 10 also shows the moderating effect of 
the stronger wind. Stronger winds increase the mixing which miaous la 
deeper layer to cool down or warm up, thus decreasing the amplitude of 
the diurnal temperature cycle. We see in Fig.ii that the amplitude of 
the diurnal variation of the roughness angle is also reduced by the 
stronger wind. The shape of both curves is similar. We observe a ra- 
pid decrease ina during the first 6 hours after the time of maximum 
temperature, However, the surface cooling rate has reached its maxie 
mum intensity near 6 PM and is decreasing rapidly between 6 PM and 8 
PM. This reduces the growth rate of % which reaches a secondary mae= 
ximum near 7 to 8 PM. & decreases after the decrease in the cooling 
rate because the level of maximum cooling rate shifts higher, desta- 
bilizing slightly the layer close to the ground. However, after a 
while a deep layer above the ground becomes so stable as to permit a 
continuous increase in % , which is observed between 11 PM and 6 AM. 
After 8 AMX exhibits a tendency to decrease, However, we observe 
a large oscillation with minimum near 9 AM and maximum near 11 AM. 

This is similar to the overnight maximum in A. The maximum warming 
rate occurs near 9 AM. Due to enhanced turbulence under unstable 


conditions, the slight decrease in the warming rate after 9 AM is 
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Figure 9. Diurnal cycle of the surface potential temperature. The 
different curves represent the followings 

— Zo 71 em, Ve = 10 m/sec, soil thermal conductivity= low. 

= 10 m/sec, soil thermal conductivity = high. 


high. 


“@— Zo= 100 cn, Ve 


S—2Z5= 1 em, Ve = 10 m/sec, soil thermal conductivity 


= 1 em, V, = 20 m/sec, soil thermal conductivity = high. 


The hour in this and the subsequent graphs igs expressed in 


Solar time related to the mountain standard time by: LST = MST~0.6, 


Sunrise occurs at 348 AM LST and sunset at ged PM LST. 
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Figure 10, Warming rate at the surface. The four curves represents 


— 2,71 cm, V, = 10 m/sec, soil thermal conductivity = low. 


g 


ree Qe 100 cm, V. = 10 m/sec, soil thermal conductivity = high. 


g 


oe) ee ue = 10 m/sec, soil thermal conductivity = high. 


ea ee OU ae Au m/sec, soil thermal conductivity = high. 
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Figure 11. Diurnal variation of the roughness angle. The various 


curves ars3 


—— 2,= 1 em, V, = 10 m/sec, soil thermal conductivity = low. 


g 
—~ % = 100 cm, V, = 10 m/sec, soil thermal conductivity = high. 


amg a ilactts ty 40 m/sec, soil thermal conductivity = high. 


“"77 Zo= 1 om, Vp = 20 m/see, soil thermal conductivity = high, 
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almost immediately reflected in “which increases suddenly. The oscil- 
lation is certainly amplified by the finite grid because the inversion 
height is not lifted continuously but by discrete steps. Therefore, 
even if there are some physical reasons for the oscillation in c& 
near 10 AM, the amplitude of the phenomenon is certainly much smaller 
than indicated in Fig. i1. 

Figs. 12 and 13 represent the diurnal cycle of the vertical 
distribution of ® for strong and moderate geostrophie winds, respectie 
vely. With moderate geostrophic winds the nocturnal cooling is only 
0.03 °C at 240 m and 0.1 °C at 166 m. However, with strong geostro- 
phic winds the overnight cooling at 166 m, 240 m and 670 mis, respec- 
tively, 12.3 °C, 6 °C and 0.1 °C. This is a reflection of the increa- 
sed turbulence and deeper mixing layer. One interesting feature of the 
model is that it predicts a delay between the time at which the surface 
temperature starts increasing after sunrise and the time at which the 
atmosphere becomes unstable in the lowest levels. Near sunrise the 
atmospheric heat flux is positive and it will take some time before 
it changes sign. The increase in surface temperature destabilizes 
the lower layer which increases the turbulence and the coefficient of 
eddy diffusivity and therefore permits a slower decrease in the heat 
flux. That delay is about one hour for the moderate wind case but 
reaches two hours in the case of strong winds. The surface tempe-~ 
rature warms up by 07 °C in the moderate-wind case and by 2.4 °C in 
the strong-wind case before the lower atmosphere becomes unstable. 

The Kansas experiments have shown an average delay of one hour bet- 
ween sunrise and the appearance of an unstable lapse rate at 6m 


(Wyngaard; 1973), This agrees very well with our results if we 
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Figure 12. Diurnal cycle of the vertical distribution of potential 


temperature for the case 2. = 4 em, V, = 20 m/sec and high soil conduce 


g 
tivity. The units are [°K]. 
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Figure 13. Diurnal cycle of the vertical distribution of potential 


temperature for the case Cra a 1 cm, Ve = 10 m/sec and high soil 


conductivity, The units are [ox]. 
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suppose that the average geostrophic wind during these experiments was 
less than 10 m/sec. 

The heat flux is represented in Figs. 14% and 15 for the cases 
of high and moderate winds. We observe that wO@ is more than doubled 
by doubling the geostrophic wind speed. The phase difference already 
noted in the surface value of wo between the two simulations is also 
apparent up to 50 m We observe also that the inversion height is 
lifted more rapidly with the stronger winds after the instability 
has been established in the lower grid points. 

By comparing Fig. 16 with Fig. 17 we can say that tho maxi- 
mum value for Ky is increased by a factor of 2 to 4 and that the height 
at which that maximum occurs is shifted upwards by about 50 m. There» 
fore, wO is larger with stronger winds, even if the temperature gra- 
dients are smaller, because of the much larger values for K, in 
that case, The general shape of both curves for K, is similar: 
very high values during near neutral and unstable situations and very 
small values everywhere overnight during stable situations. We observe 
also a marked difference between the first day and the second day due 
to the fact that the upper levels have remained stable during the 
24ehour cycle which caused a decrease in K, above 500 m. The model 
hes a tendency to produce negative values for K, under very stable 
conditions. We reset Ky Seeeeiealie to a very small number of the 
order of 10° em“/sec whenever that happens. We need that restric- 
tion in order to insure the stability of the finite-difference scheme. 
Otherwise we would get the formation of many artificial inversions 
which could easily prevent the convergence of the model. This 


problem is more serious in the moderate wind case as 5 to 6 levels 
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Figure 14. Diurnal cycle of the vertical distribution of w6 for 
the case z, = 1 em, V, = 20 m/sec and high soil conductivity. The 


units are [em ec sec"!], 
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Figure 15. Diurnal cycle of the vertical distribution of wo for 


the case z 


pr ceecs ¥o.= 10 m/see and high soil conductivity. The 


units are|om Ae sec! |. 
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Figure 16, Diurnal cycle of the vertical distribution of the coef- 
ficient of diffusivity for temperature for the case z, = 1 cm, Vg s 


10 m/see and high soil conductivity. The units are [10* en@/secl, 
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20 m/sec and high soil conductivity. The units are [10% en/sec Jo 
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are computed with negative values by early morning. With strong geoe 
strophic winds only one or two levels exhibit that tendency. We feel 
that allowing negatives values for K, would not change significantly 
the results while compromising the convergence of the iterative process. 
This is the justification for not allowing K, to become negative. 

Figs. 18 and 19 show the diurnal variation of U for moderate 
and strong winds. Both graphs are similar in shape. The maximum wind 
gradient is near the ground under unstable or neareneutral lapse rates, 
The wind gradient decreases overnight in both the upper boundary layer 
and the lower part of the surface layer. The region of large overnight 
wind gradient is concentrated between 20 m and 50 m for the moderate 
wind case and between 20 m and 166 m for the strongewind case, Between 
about 100 m and 300 m, U has maxima near 12 PM and 12 AM and minima near 
5 AM and 5 PM. The overnight maximum is caused by an inertial oscilla-= 
tion which takes place some time after the initial coaling when the 
friction force decreases above 50 m within a short period of time, 
from a relatively large daytime value to an extremely small value. 
There is a tendency for the wind to "overshoot" its geostrophie value 
which is caused by the non-zero time derivative term in the equations 
of motion. Fig. 20 illustrates the diurnal cycle of the total wind 
speed at different heights for the case of strong geostrophic winds. 
Near the ground the first minimum occurs near 8 PM, just after the time 
of maximum cooling rate. The decreasing cooling rate after 6 PM permits 
a small recovery in the surface wind speed which reaches a secondary 
maximum near 12 AM. The stabilization of the cooling rate after 10 
PM causes a further decrease in the surface wind until sunrise. After 


sunrise, the surface wind increases to reach its maximum near the time 
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Figure 18. Diurnal cycle of the vertical distribution of U for the 


case of z, = 1 em, Vz = 10 m/see and high soil conductivity. The 


units are [m/sec] ° 
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Figure 20. Diurnal eyele of the total mean wind for the case Zy = 1 cmy 


= 20 m/see and high soil conductivity. The various curves represent 


the following heights: 
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of maximum temperature. With strong winds the maximum amplitude of the 
double-maximum-wind eyele is reached near 404m. The overnight maximum 
near 12 AM is due to the inertial oscillation and the maximum near 12 
PM is due to the influence of surface heating, 

Figse 21 and 22 represent the diurnal cycle of V. In the sur- 
face layer, V increases generally overnight causing an increase in the 
roughness angle. V reaches a secondary maximum near 8 PM in the vicie 
nity of 40 m for moderate winds and of 80 m for strong winds. V shows 
also a daytime maximum near 11 AM. The region of overnight maximun 
gradient for V is slighly higher than the region of maximum gradient 
for U. Above that region, V undergoes a double cycle with maxima near 
4% PM and 10 AM and negative minima near 2 AM and 2 PM. The influence 
of the higher wind speed is to spread out the velocity gradient over 
a larger region and to shift the position of the maximum winds higher 
in the atmosphere, 

Figs. 23 and 24% represent the diurnal cycle of e. Wea have 
a rough estimate of the constant-stress layer by looking at the layer 
in which e does not vary by more than 20%. For the moderate-wind case, 
the height of the constantestress layer is near 120 m under neutral 
conditions and could be as high as 400 m during unstable situations. 

It seems more difficult to define a constantestress layer overnight 
because the stress reaches its maximum value above the ground, near 7 me 
In any case, it seems that the lowest overnight constant-stress layer 
was around 20 my and happened just after sunrise. The strong wind 

ease is similar with the turbulent energy in the surface layer being 
twice as large during daytime and 23 times during nighttime. The height 


of the constant-stress layer varies from its value of 240 m 
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Figure 21. Diurnal cycle of the vertical profile of V for the case 


2, = 1 em, V, = 10 m/see and high conductivity. The units are 


[m/sec]. 
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Figure 22. Diurnal variation of the vertical distribution of V for 


the case 2, = 1 em, V, = 20 m/see and high soil conductivity. The 


g 
units are[ m/sec). 
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Figure 23. Diurnal variation of the vertical profile of e for the 


case z, = 1 cm, V, = 10 m/sec and high soil conductivity. The 


units are [om*/sec?]. 
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Figure 24. Diurnal cycle of the vertical distribution of e for the 
case 2, = 1 cms Vy = 20 m/sec and high soil conductivity. The units 


are [em?/sec?]. 
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under neutral conditions to above 320 m under unstable lapse 
rate and goes down to 120 m overnight. The maximum daytime value of - 
the constant-stress layer depth depends strongly on the depth of the 


unstable layer. 


6.2.2 Effect of the Roughness Height. 
We go back to Fig. 9 and note that the temperatures over 


rough terrain always remain lower than the temperatures over smooth 
terrain. The initial temperature difference is 1.4 °C, it anovantss 

to 1.8 °c by sunrise and to 2.7 °C by the time of maximum temperature. 
The greater cooling over the rough terrain during the od enttine is 
contrary to theoretical expectations because we expect that with 

the enhanced turbulence over the rough terrain, a deeper layer will 
cool down, thus increasing relatively the surface temperature. Fig. 

10 shows that the cooling rates over smooth and rough terrain are 
elneee identical until 6 AM. A more detailed analysis of the data 
reveals that between 2 PM and 8 PM, the air over the rough terrain 

cools dow more slowly in agreement with the theory. However, this 
situation is reversed between 8 PM and 4 AM. When we compare Fig. 25 
with Fig. 13, we observe that the temperature gradient in both cases 
accumulates between 20 and 100m. For z, = 1 em, the temperature diffe~ 
rence near sunrise between 20 m and 50 mis 10.6 ccs and between 50 m 
and 100 mis 7 °C. Similarly, for ona 100 em, the difference is 6.5 °C 
between 20 m and 50 my and 9.7 °c between 50 m and 100 m. Therefore, 

a deeper layer is cooling over the rough terrain as indicated by the 
spreading out of the temperature gradient. Now, we compare the heat 


fluxes in the atmosphere as represented in Figs. 15 and 262 In the 


oadal sidatan: cetat « 08e arectn, ad 2th 
*o suluy. omtiq¢ed mmt kon ol Sight ® ts of a 
ait bh Wrqeb ait megs gett one wwe 


t 


‘aemng eaveettse tend etd Spee ston ban @ et tied 99. a 
ee AS gail ewe adsei emrl | 
Labo d ae (oP ne L2 eoneTeN LAD exmhemmgret ip Mind ef? 
otenttee? msmkea So web dade tn 
} et ent ifdete ate galrak sisi Gace a 

tty Ana eR: ey ‘vue anndtesveqne Leottexood} : 
shige Habel hae <a ot var: man Sea ope 

alt geri aagease « eagelai ts ehebibaten'gtouwnon sty tt i 

| wa nase eH ain, ter yrtien wide Saat 5 in 2 
ial itt vo aac ttle tom A. si nic aetna 
cep gps 0 ih se nt 8 eS yt 

"geet ROTOR sehr aah aac Pohemertan et Side eet tna le 200 
ee ee 10, + bas BBD creditor sei at a 

ama tied 2.200 Maran eeoet Wh dni? avewede. ae oft sit athe’ 

ee a ee a0 Stage sie 
wi ef a? PT 8 M08 am» 28 stielicd qs tiles tanita 
Paihia hae yn POE# gee wp? aybradiett iy ohh ter be ; 
ou reed at Ot busy qaimmndeet °° S.0 tes lO Bee OS oeeeded ~ 
a4 & betae et 2A «ta et Ad 197 “nilpoe by saya Pare * 
ee ee eit te sue yabhawagh 4 


bad 


haati “48 Ppa: ¥ 


bees 
A ai!) temp < Ab Wetter 24 — nit nz ‘pout? y 
~ _ Ye) ae iy 
: ~ vi yf ; " @ P 


122 


z (m) 
1832 
1655 
1478 
1302 
1127 
954% 
784 
617 
455 
301 
163 
Bye 


8.8 


EE iG aR oy “Aa a” TRO SCerige jor ae 2 
PM AM PM 


Figure 25. Diurnal cycle of the atmospheric potential temperature 


for the case z, = 100 em, Ve = 10 m/sec and high soil conductivity. 


The units are [ x}. 
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Figure 26. Diurnal cycle of the vertical profile of wO@ for the case 


z= 100 em, V_ = 10 m/sec and high soil conductivity. The units are 
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lower boundary layer, the heat flux is about twice as large over rough 
terrain between 4 PM and 8 PM. After 8 PM, the difference between the 
two heat fluxes diminishes to become negligible by sunrise. The influ- 
ence of the larger roughness height is reflected in the fact that the 
heat flux has a larger value over a deeper layer over rough terrain 
than over the smooth terrain. There is also more diffusion over rougher 
terrain as demonstrated by the higher values of Ky in Fig. 27 as come 
pared to Fig. 16. 

Comparing Figs. 28 and 23, we note that the turbulent nares 
in the surface layer is increased by 50% over rough terrain during 
daytime, and only by 25% overnight. 

Figs, 29 and 18 show that the diurnal cycle of U is similar 
over smooth and rough terrains. The overnight wind gradient which was 
concentrated between 20 m and 50 m over smooth terrain is now spread 
between 20 m and 100 m over the rougher terrain. The double-wind- 
maxima feature above 100 mis also very evident over rough terrain. 

The diurnal cycle of V over rough terrain is very similiar to the one 
over smooth terrain as demonstrated by the comparison of Fig. 30 with 
Fig. 21. The rough terrain induces more turbulence which causes a 
spreading of the wind gradient and a shifting of the maximum wind speed 
higher. The only other major difference is in the roughness angle 
which is illustrated in Fige 11. The roughness angle is generally 

10 to 15 degrees greater over the rougher terrain. The oscillation in 
XH after sunrise has a larger amplitude over the rough terrain and 
reaches 40°, 

Therefore, the rougher terrain creates more turbulence at 


all times even if the difference in turbulent energy and heat flux 
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Figure 27. Diurnal cycle of the vertical profile of K, for the case 


2, = 100 em, Ve = 10 m/sec and high soil conductivity. The units are 


[10* em/sec]. 
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Figure 28, Diurnal cycle of the vertical profile of e for the case 


Zo = 100-om,; V-.= 10 m/sec and high soil conductivity. The units are 


[em?/seo], 
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Figure 29. Diurnal variation of the vertical profile of U for the case 


Zz. = 100 cm, V, = 10 m/see and high soil conductivity. The units are 


[m/sec]. 
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Figure 30. Diurnal variation of the vertical profile of V for the case 


Zz, = 100 cm, V_..= 10 m/see and high soil conductivity. The units are 


[m/sec]. 
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is reduced overnight due to the stability. The greater cooling over 
the rougher terrain between 8 AM and 4% PM remains unexplained. A plau- 
sible explanation would be that the two simulations have not converged 
sufficiently to their true values at each time step, creating an arti- 
ficial difference between the two simulations which has persisted overe 
night. Therefore, it seems necessary to require more accuracy than 10 om” 
/sec overnight because the maximum value of Ky is reduced by an order 


of magnitude overnight which reduces considerably the relative accu- 


racy with which we solve the equations. 


6.2.3 Effect of the Soil Conductivity. 
Fig. 9 shows that the effect of decreasing the soil thermal 


conductivity is to amplify the diurnal temperature wave at the sure 
face of the earth. The maximum surface temperature is increased by 
18.5 °C and the minimum by 2.8 °C. As indicated in Fig. 10, the ware 
ming rate is twice as large over the poorly-conductive soil, The graphs 
of Hand wo over the poorly=conductive soil are Figs. 31 and 32. These 
are very similar to Figs. 13 and 15 which represent the same va- 
riables over a good conductive soil. The only major difference is in 
the intensity of the cooling which is more pronounced over a soil with 
low thermal conductivity. Fig. 33 reveals the influence of the ampli- 
fication of the diurnal temperature wave on the atmospheric diffusivi-g 
ty. Over a poorly conductive soil, Ky 4s larger under unstable cone 
ditions and smaller under stable conditions. 

The behavior of e over a soil with low conductivity is repre- 
sented in Fig. 34 and is basically the same as over a soil with high 


conductivity as represented in Fige 23. Due to strong cooling, the 
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Figure 31. Diurnal cycle of the vertical profile of the potential 
temperature for the case 2, = 1 cms Ve = 10 m/sec and low conducti- 


vity in the soil. The units are [°K]. 
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Figure 32. Diurnal cycle of the vertical profile of w © for the case 


z, = 1 om, V, = 10 m/sec and low thermal soil conductivity. The units 


° g 


are [em-°C/sec). 
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Figure 33. Diurnal cycle of the vertical profile of Ky for the case 


ona 1 em, v, = 20 m/sec and high soil conductivity. The units are 


[cm*/sec]. 
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Figure 34, Diurnal cycle of the vertical profile of e for the case 


a 1 em, ye = 10 m/see and low thermal soil conductivity. The units 


are [cm?/sec* |. 
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Figure 35. Diurnal cycle of the vertical profile of U for the case 


Dear 1 em, V 


[m/sec]. 
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Figure 36. Diurnal variation of the vertical profile of V for the 


case 2, = 1 cm, v, = 10 m/see and low soil conductivity. The units 


are [m/sec]. 
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minimum near 6 PM is 15% lower over a poorly-conductive soil. However, 
@ recovers after 12 AM and has almost the same profile as over soil 
with a high conductivity. 

Figs. 35 and 36 represent the U and V profiles over a poorly. 
conductive soil, These are similar to the profiles over a good-conduc- 
tive soil as represented in Figs. 18 and 21. The overnight minimum in 
U within the surface layer is reduced by about 20% over a poorly.conduc- 
tive soil. Above the region of maximum wind gradient, the double-wind- 
maxima cycle is still observed but with an amplified amplitude, As 
shown in Fig. 11, the variations inaare amplified over a poorly .conduc- 


tive soil, with higher maxima and lower minima, 


6.2.4 Examination of the Other Variables. 

We will examine very briefly the behavior of the other 
atmospheric variables. Only one specific case will be studied: ors 
1 cm, vn = 20 m/sec and high soil conductivity. This case was chosen 
for the following reasonss 
a) the profiles of most of the variables are smoother with strong 
geostrophic winds. 

b) The overnight gradients in @®, U, V and e are spread out over a 
few more grid points. 
¢) K, becomes negative at fewer grid points. 

The diurnal variation of is fairly regular as shown in 
Fig. 37. The diurnal variation is negligible close to 2,, reaches 8% 
of the mean value at 21 m and exceeds 50% of the mean value near tho 


top of the boundary layer. Due to its definition, | is always increae 
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Figure 37. Diurnal variation of the vertical profile ot J for the 


case Zo 


A 


= 1cm, V 


g 


= 20 m/sec and high soil conductivity. 


The units 


‘ sit A tet, 


AR 05) horses e 


ay — yo 


et er ik iiss Aa a Saar bE: ™t ie aa or 
a aot Pt. _ ‘ ia wos 
vay 


“ 7 ‘ se sh) my 


| wie . separ ome rt Oi vind sesh Sh 


inti 


138 
sing upwards. The delay of 6 to 8 hours between the time of change in 
stability in the lower boundary layer and the time of minimum for / 
indicates clearly that the formula which computes f depends to a large 
degree on the value of e in the uppermost levels. 

As illustrated in Fig. 38, O2 reaches an overnight maximum 
in the lower boundary layer just after the time of maximum cooling 
rate and decreases until the air becomes unstable near 6 AM. Between 
8 AM and 12 PM there is a rapid increase in 92 which corresponds to the 
increase in the warming rate. Under unstable conditions the es 
value of Q2 is always near the surface. However, the temperature gra~- 
dient accumulates overnight between 20 m and 100 m. This is reflected 
in the large overnight maximun moe in that region, 

The diurnal cycles of u6 and V6 are show in Figs. 39 and 
40, Equations (2.195) and (2.196) indicate that these variables are 
proportional to the product of the appropriate wind eraatens by the 
temperature gradient. Therefore, under stable conditions, these vae 
riables are positive if the wind increases with height and negative if 
the wind decreases with height. As the wind and temperature gradients 
accumulate overnight between 20 m and 100 m, the maximum value is 
reached in that region. 

Figs. 41 and 42 illustrate the stress in the x and y direc- 
tions. In the lower boundary layer, uw has a large negative values 
during daytime and a smaller negative value during nighttime. This 
is another illustration of the effect of the decreasing turbulence 
overnight. Vv w decreases between 2 PM and 6 PM but increases after 
6 PM. Therefore, the overnight increase in roughness angle is accom~ 


panied by an increase in the stress in the y direction. These vae 
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Figure 38, Diurnal variation of the vertical profile of Q° for 


the case z, = 1 cm, Vp = 20 m/sec and high soil conductivity. 


The units are [ oc?]. 
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Figure 39, Diurnal cycle of the vertical profile of u@ for the case 


z= 1m, V, = 20 m/sec and high soil conductivity. The units are 
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Figure 40, Diurnal cycle of the vertical profile of VO for the case 


Zz 


9° 


1 cm, Ve 


Lom-°C/sec)s 


= 20 m/sec and high soil conductivity. Tho units are 
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Figure 41. Diurnal variation of the vertical profile of u w for the 


case 2) = 1 em, Me = 20 m/sec and high soil conductivity. The units 


are [ m?/sec?]. 
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Figure 42, Diurnal cycle of the vertical profile of v w for the case 


z= 1 om, V, = 20 m/see and high soil conductivity. The units are 


[m?/seo?]. 
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riables depend on the vertical gradient of the mean wind components, 
Fig. 43 illustrates the behavior of K,, which has lower daytime 
values than Kye However, Ke has a smaller tendency than Ky to develop 
negative values overnight. 
Figs, 44% to 46 illustrate the behavior of ie, me and 3 whose 
sun is simply oe”. Therefore, their behavior is very similar to the 


behavior of e. They reach a maximum during daytime and a mini- 


2 and Z 


mum overnight. The approximate daytime contributions of uv, Vv 
to eo are, respectively, 45%, 32% and 23%. Fig. 47 represents the 
behavior of uv. Near the surface, the minimum is reached just after 
the time of maximum cooling rate and the maximum is obtained just af- 
ter the time of maximum warming rate. This variable is proportioc- 
nal to the product of the gradients of U and V which explains the poe 
sitive sign in the lowest 50 m and the negative sign in the region 


above 100 m. 


6.3 Comparison of the Soil Temperatures. 
Figs. 48 to 51 represent the diurnal variation of the soil 


temperature for the four simulations, The features common to all graphs 
ares 

j- There is a time lag between the upper and lower levels in the time 
at which they repectively reach their maximum and minimun temperatures. 
2 The temperature variation below 0.5 m is very small, This indicates 
that the soil depth of 1 mis appropriate for the numerical simulations 
for the periods of the order of one days 

3. The temperatures near the bottom increase almost continuously 


because the surface minimum temperature is too warm. The other levels 
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Figure 43, Diurnal cycle of the vertical variation of K, for the 


case z. = 1 em, V, = 20 m/sec and high soil conductivity. The units 


° g 
are [10° en’/sec}. 
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Figure 44, Diurnal variation of the vertical profile of u% for the 
caso z. = 1 omy bes = 20 m/sec and high soil conductivity. The units 


o 
are [m2/sec}. 
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Figure 45, Diurnal variation of the vertical profile of v* for 


the case z, = 1 em, V, = 20 m/sec and high soil conductivity. The 


units are [m®/sec* ]. 
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Figure 46, Diurnal variation of the vertical profile of w* for the 
case Z = 1 cm, Vy = 20 m/sec and high soil conductivity. The units 


are [ m2/sec?] : 
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Figure 47. Diurnal variation of the vertical profile of uv for the 
case 2, = 1 em, V, = 20 m/sec and high soil conductivity. The units 


are [m?/sec?| ; 
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Figure 48, Diurnal variation of the soil temperature profile for the 


case %, = 1 em, Vy = 10 m/sec and high soil conductivity. Tho units 


are (°x]. 
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Figure 49. Diurnal variation of the soil temperature profile for the 


case Zz) = 1 em, Me = 10 m/sec and low soil conductivity. The units 


are (° |, 
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Figure 50. Diurnal variation of the soil temperature profile for the 


case 2 = 100 em, V, = 10 m/see and high soil conductivity. The wits 


are [°K] ° 
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Figure 51¢ Diurnal variation of the soil temperature profile for the 


case 2, = 1 cm, ve = 20 m/see and high soil conductivity. The units 


are [ °K |. 
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are not perfectly cyclical. The second-day maximum is higher than that 
of the first day due to the decreased soil heat flux from one day to 
the other. The surface temperature could be lowered by modifying 
either the transmissivity factor or the albedo. In this way we could 
obtain more cyelical diurnal variations within the soil layer. 

We now compare a poorly-conductive soil with a good-conduce 
tive one, The poorly-conductive soil has: 
ij a larger temperature variation in the upper layers. 

2- A smaller temperature variation in the lower layers, 
3- The time lag between the maximum or minimum temperature in the 
upper levels and the lower levels is considerably increased, 

All these differences are expected because we expect that 
the temperature wave will not penetrate a poorly-conductive soil 
as deeply as it would a good-conduective soil. This tends to cone 
centrate the temperature gradients very close to the surface in the 
case of a poorly-conductive soil. 

There is almost no difference between the soil temperature 
profiles due to a change in roughness height. The temperatures are 
slightly warmer in the case of 2, = 1 em. A detailed calculation of 
the soil heat flux during the period 8 PM to 2 AM shows that the heat 
flux from the soil is about 10% higher with 2, = 100 em. Therefore, 
the atmospheric heat flux and the soil heat flux are both larger 
overnight in the case of Beis 100 em, which should give a warmer sure 
face temperature in this case. The number of itera- 
tions taken by the model with 2, = 100 em was approximatively 2 to 3 
times the number of iterations required by the model with 2, = 14 em, 


This seems to have caused a significant difference in the convergence 
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Table 5. Comparison of the surface potential temperature between 


the four simulations of the urban heat island. 


model # 
parameter 


minimum G)(2,) 289.5911 292.5137 - 287.9182 -292.5120 
maxinun Gz) 311-9797 310.2654 309.2634 330.2451 
time of minimum yO ay ea ae NT yD py 


fame (of) maxi mum 220 py 220 py 220 py 


minimum wO (z,) 
maximum w@ (25) 
time of minimun 


time of maximum 


time of wQ (z,) = 0 


Surface warming rates 
minimun 
maximum 

time of minimum 


time of maximum 


Zo Lem] 
v, (n/sec ] 


soil conductivity 
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of the surface temperature for both simulations. The model usually 
converged without any oscillation when we used OMEGA = 0.7, which 
increases the probability of a build-up of systematic errors as i3 


the case here. 


6.4 Assessment of the Urban Heat Island Effect. 

We compare the quantities related to the behavior of the sure 
face temperature for the four simulations in Table 5, Usually 
the city has a larger roughness height than the surrounding peuiteysiien 
The city has also a higher soil thermal conductivity because of the 
abundance of concrete in the city and also because of the fact that 
clay, a poor thermal conductor, is often present in the countryside. 
Therefore, we expect a smaller diurnal amplitude over the city because 
of the increases in roughness height and soil thermal conductivity. 
The results seem to indicate that the countryside would be warmer 
than the city during daytime. However, another major factor coune 
teracts the strong effect of the soil thermal conductivity: evapora- 
tion. There is generally a very important decline in evaporation 
over the city. This permits a large part of the solar energy to 
go into latent heat over the countryside. The difference in evapora- 
tion between the city and the countryside needed in order to obtain 
a warmer city depends mainly on the differences in the soil 
properties. 

When the effects of evaporation and soil thermal properties 
nearly balance each other, the influence of other factors may become 
important in a specific case when we wish to determine which of the 


city or the countryside will become warmer during the daytime. The 99 cm 
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increase in roughness height cools down the city by 2 °C during day- 


time. The influence of an increase in the geostrophic windis to reduce 
the amplitude of the diurnal temperature cycle and consequently to 
reduce the urban heat island effect in absolute valus. 

There is no significant difference in the times of maximum 
and minimum temperatures for all four simulations. However, the time 
of change in stability in the lower atmosphere is retarded by strong 


winds and by high soil conductivity. 
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CONCLUSION 


A numerical model using a more sophisticated approach to ture 
bulence modelling than the ordinary K-theory has been applied with some 
success to the problem of the urban heat island. The unidimensional 
model basically describes the interaction between a soil layer and 
an atmospheric layer. Three factors were examined: roughness height, 
soil conductivity and wind speed. The model has demonstrated that the 
soil conductivity was an extremely important factor in determining the 
amplitude of the diurnal temperature cycle. The maximum temperature was 
considerably increased over a poorly conductive soil, which is 
often found in the countryside. Increased roughness height 
also reduced the amplitude of the diurnal temperature wave by a small 
amount. When we associate a high roughness height and a high soil cone 
ductivity with a city, we are led toe the conclusion that the city is 
considerably cooler than the countryside during the daytime. The cone 
terbalancing effect which was not studied is the difference in evaporae 
tion between the city and the countryside. If we suppose a dry city and 
a moist countryside, we expect that the daytime rural temperature will 
be cooler. 

Many problems are associated with the numerical solution of 
the equations. An implicit finite-difference scheme was used in order 
to provide numerical stability and to permit a time step of 10 minutes. 
The equations have to be written very carefully in finite-difference 
form in order to reduce the round-off errors to which the temperature 


computations are very sensitive in the upper boundary layer. Double 
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precision had to be used for the mean quantities in order to 
obtain reasonable profiles above 1000 m. In some cases systematic 
errors have built up due to insufficient accuracy in the itera- 
tive procedure within each time step. The accuracy which we need in 
order to obtain good results must be determined by trial and error. 
The disadvantage of requiring more accuracy is the increase in compu- 
ting time. 

The most obvious improvement to the model is to include the 
moisture computations which present a few supplementary problems. We 
nesd a balance equation between the soil and atmospheric moisture. 
There is also the problem of water vapor build-up in the atmosphere, 
The main source of water vapor in the boundary layer is the evaporae 
tion at the surface which is enhanced during daytime. The water vapor 
is also transported fairly high in the atmosphere under unstable lapse 
rates, However, the main sink of water vapor is not the overnight 
condensation at the surface which involve only a shallow layer close 
to the ground under stable conditions. The main sink is the formation 
of clouds and precipitation. Therefore, we would need to incorporate 
some cloud physics into the model or we could remove artificially 
any excess water. If wa include condensation in the atmosphere, we have 
to introduce a source term in the equation of temperature due to the 
release of latent heat. This would imply also a modification of the 
solar energy received at the surface of the earth due to reflection by 
the clouds and to Sas Fea within the cloud. The infrared balance 


would also be affected. 


Other improvements may include an albedo model which permits 
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a diurnal variation of the albedo with the solar angle. Radiative flux 
divergence could be important in the lowest 100 m over a city. The dise 
advantage of the roughness height concept over a city is that the value 
associated with very dense and tall buildings is of the order of a few 
meters. Therefore, we do not model the lowest meters over the city 
- although often this is the layer in which we are interested. We could use 
twoe- and three-dimensional models in which some care has to exercised 
in order to model properly the change in roughness height. As pointed 
out earlier, the forced-convection parameterization of some of as 
turbulent quantities gives a different result than the one predicted 
in the free-convection limit. Therefore, we must reconcile the two 
points of view. We would also like to model more complex soil layers in 
which the amount of water influences greatly the soil conductivity, 
Also, we could investigate the effects of two or more types of soil 


layers superposed one upon another, 
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APPENDIX A 
LOG-LINEAR ATMOSPHERIC GRID 


As explained in Chapter III the best transformation of coor- 


dinates in the lower boundary layer is, 


ys act In|(2#29)/29] /kotz/ fo} (A.1) 


= log-linear coordinate 


ac 


where, 


ordinary vertical coordinate (cm) with origin at Zo 


= roughness height 


nN 
| 


= von Karman's constant 


=> 
°o 
it 


length fixing the maximum size of the eddies 


constant chosen for our convenience. 


AC 
The grid is defined by the following statements, 
- The maximum number of ¥ levels is 40 


- Y and z are zero at Zo 


Ly 


the top of the boundary layer is at 3000 m and correspond to Y= 39 
- the number of grid-points is fixed at 40 and the index k represen- 
ting the height of the grid-point is related to k= Y+ 1. 


The constant AC is determined when we assign a value to Zo 


and He 


ii 
‘c xe ( 


AC = 7 ———sonntoFig + 300000 A?) 
eae aetecs ats 


° 


During the transformation of coordinates various derivatives 
of Y with respect to z appear; these derivatives can be derived exac- 


tly from equation (A.1). 
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Therefore we only have to solve the equation for z at each 

Y level, There are three iterative ways of solving (A.1) for 2. 
Method 1. Near 2, the grid is almost logarithmic because of the small 
contribution from the linear term in that region. The most rapidly 
converging formula is obtained by isolating the term inside the lo- 
garithm, 

Ze 25¢exp[Kg( Y /ACe2/  ,)]- 1} (A.5) 
Method 2. Near the top of the boundary layer the grid is almost linear 
because the logarithmic term has a small influence there. Therefore 
we isolate the linear term, 

z = fo} ¢/ac-In[(2#29)/24] ig} (A.6) 
Method 3. If formulas (A.5) and (A.6) do not give reasonable re- 
sults we can use formula (A.1) and do the iteration on Y instead of 
on z In this case we readjust the value of z until we get close 


enough to the integral value for a given level. 


A.1i Doatermination of IKM. 

Firstly we replace ¥ by k =4+ 1 because most computing 
systems do not accept an index 0. The level IKM is the level repre- 
senting the upper limit of the applicability of (A.5)« Close to the 
roughness height the initial guess will be purely exponential, 

2, = nof exp[(Keg=1)/ac]-t } (A.7) 
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The initial guess is always greater than the true value for 
z because we have neglected a negative number in the exponential, name- 
ly -2;// 5+ At some level the initial value of z will be large enough 
to cause the term inside the exponential in (A.5) to become negati- 
ve: that level can be considered as the first one where the exponential 
initialization is inappropriate. Therefore, KM is the level immediate- 
ly below the first level where, 

(ret) {icra 2/1 G0 (A.8) 


A.2 Initialization of the Other Levels. 
We use the assumption of a linear spacing to initialize all 
the other levels, 
2; = YxXtop of the B.L.)/(maximum ¥) = (k-1)x300000/39 (A.9) 

That initial guess is very good for the top levels where 
the spacing is indeed almost linear but there is a critical region 
where neither the exponential spacing nor the linear spacing are 
expected to be good. One way of determining that region is to input 
z, into (A.6) in order to get a second guess 2’, 

2! “25 4 pac - In] (24429)/ 29] [oh (A.10) 

The overestimation in Zs in the lower levels can cause 2! 
to become negative which would stop the iterative procedure because 
the next guess would require the estimation of a logarithm of a ne- 
gative number, Therefore, each time that z' is negative we readjust 
zy toa smaller number 

Z, = 2, /B (AS) 


where B is a positive constant greater than 1. 
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A.3 Accelerated Gauss-Seidel Process, 

We use the accelerated Gauss-Seidel process which consists 
in refining our initial guess by giving a weight both to the last 
approximation and to the new computed value. For example, we call, 

Z = last value of 2 
TEMP = computed value for z 


TEMP. Z = difference between two successive iterative values 


TP 
of z 
Z2* = Z + TP OMHGA = new approximation to z 
OMEGA is an acceleration factor which gives a weight both to Z and 
TEMP, It can vary from 0 to 2 and, 
OM2GA = 1 implies that the computed value of z is TEMP. 
OMBGA 5 1 implies acceleration, appropriate for the cases in which 
Z* approaches a value steadily without oscillations. 
OMBGA < 1 implies deceleration, appropriate for the cases where 
Z* oscillates around the real value of 2. 

Unfortunately it does not seem possible to devise a general 
formula for the optimum value of OMHGA in our problem. Some of the 
difficulties are, 

(a) close to 2,, OMEGA = 1 is fairly appropriate because the oscilla- 
tions around the real value are damped very rapidly. 

(bo) Near IKM large oscillations prevail due to the fact that the ini-e 
tial guess is poor in that region and also that neither (A.5) nor (A.4) 
are really appropriate thore,. 

(¢) The level IKM is dependent on the roughness height 2, and to a 


lesser extent on). 


(d) Above IKM the iterative formula changes and most of the diffi- 
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culties arise near IKM where OMHGA < 1 is needed for convergence. 
Near the top of the boundary layer the convergence is rather rapid with 


OMBGA=1. 


A.4Y Detailed Analysis of OMHGA in a Specific Case. 
With OMEGA = 1 two successive iterations will always give 


one value above and the other one below the actual value of 2. This 
phenomenon is due to the form of the equations. When the formulas are 
applied in their appropriate regions the oscillations are damped very 
rapidly. It would be difficult to determine the best OMEGA a priori 
for all the grid points. However, after the solution has been obtai- 
ned we can compute what the best OMHGA should have been. The best 
OMSGA is defined by 

OMBGAt = (Z'=Z)/(TEMP=Z) 
where, Z* = actual height of the grid point 

Z = last iterative value for z. 

The best OMEGA is therefore defined as the one which would 
give the true value in only one iteration. We will study OM@GA in the 
special case where Z, = 1 cm and J, = 2700 cm. In this case we find 
that AC = 0.2650463308, In Tables 6 and 8 we list the computed values 
of the grid-point heights at various stages of iteration. In Tables 7 


and 9 we have tabulated the computed best OMEGA at all stages of ite- 


ration. 


A.4.1 Analysis of the Iterative Process Below IKM,. 
We note that the initialization is extremely good close to 


the ground and worsens considerably near IKM where the initial value 
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Table 6 e 


stages of iteration. 
the initial values for which fomula (A.7) was used, 
cessive values are computed from (A.5). 


terms of k =Y+ 1, 


2(2) 


2 «745382866 
2 744050185 
2 0744050832 
2744050831 


z(3) 


1302789281 
13200422248 
13.00426545 
1300426538 


z(7) 


2760.433670 
1929 .065229 
2148, 690295 
2088.351785 
21044758060 
2101.284420 
2101.230814 
2101.245422 
21010241444 
2101242526 
2101.242311 
2101.242289 
2101.242295 
2101.242293 


(4) 


51.53982937 

118997609 
51219234303 
51619232702 
§1.19232713 


z(8) 


10337-88097 
27050932184 
7279 6085837 
4023201229 
6136303980 
46652759005 
5645811378 
49724118966 
5425 .952268 
51152893194 
53256744693 
5182.995555 
5279 «748775 
5213.794914 
5258.570325 
5228.131046 
5248,805082 


5240,437109 


z(5) 


195.78177570 
190.08504693 
190209731477 
190.09700948 
190.09701708 
190.09701689 


2(9) 


3872306763 
25568352926 
3725892336 
300.3897878 
37243419405 
308.9370814 
37201295099 
310.5985422 
3719393929 
310.9235327 
37192237819 
310.9854507 
37192.07383 
3109977602 
37192201448 
311.0001605 
37192.00291 


31120007420 
3719200009 


The levels are in 


z(6) 


736 60339460 
662 69963432 
674.8039290 
674.2916315 
674.3364784 
674. 3325524 
674,.3328961 
674.3328660 
6743328686 
6743328684 
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Computed values for the grid-point heights z2(k) at various 
The first entries in each columns are 
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Best OMSGA corresponding to the computed values of the EER” 


points at various stages of iteration. 


2(2) 


0.99390 
0.99845 
0.999 


2(3) 


0.9981 
0.9983 


2(6) 


0.792898 
0.783962 
0.78635 
0.785745 
0.78819 
0.7737 
0,83 
0.72 


2(4) 


0.99327 
0299328 


2(7) 


0.6679 
0.5542 
0.6261 
0.5760 
0.6092 
0.58637 
0.60172 
0259122 
0.59822 
0.5934 
0.5967 


2(5) 2(6) 
0.9978 0.9198 
0.9957 0.91948 
0.9957 0.91951 

0.91991 

2(9) 

0.2695 

0.7305 

0.2695 

007305 
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Computed values of the grid-point heights z(k) at various _ 
We have used formula (A.10) for the 


initial values (first entry in each colum) and formula (A.6) 
for the other values. 


z(10) 


75000.0000 
5087 .24000 
25843-0580 
13306.2064 
18426.7500 
15915.3529 
17045.5800 
16516.3606 
16759.6509 
16646,8531 
16698.9450 
16674, 8443 
16685.9853 
16680.0833 


16682.5 


Table 9, 


2(11) 


82500.0000 
14538.9041 
27930 «2370 
228939750 
24427 67900 
23927 6 5590 
24.087 1600 
24035.8800 
24052-3200 
24047 .0400 
24048.7400 
24048 .1970 
24048,3718 
24048 .3160 
24048.3339 
244048 , 3281 
24048 . 3303 


z(12) 


90000-0000 
2405465795 
34233 2230 


32150-3100 
31995 4144 
320326700 
3202526900 
3202528500 
32025-3300 
3202564599 
3202504291 


z( 13) 


97500.0000 
33624.0000 
41836 .5800 
40150.8100 
40468, 0000 
4OL07 «3700 
40418.9500 
4041627400 
404171640 


at some of the iteration steps. 


z(10) 


0.8340 
0.5586 
0.7306 
0.6590 
0.6945 
0.6860 


2(11) 


0.7100 
0.7700 
0.7526 
007585 


2( 12) 


0.879 
0.783 
0.811 


z( 13) 


0.8930 
0.8270 
0.8421 
0.8419 


2(14) 


105000.000 
43239 .2259 
50083 23298 
489498161 
49098.6320 
4910249959 
494102.3103 
4910224180 
49102-4011 
49402.4037 
4910224033 
49102.4034 


z( 14) 
0.865 


z(15) 


112500.000 
528938983 
5871505240 
579100400 
5801626000 
58002 4190 
580043050 
58004. 0544 
580040877 
580040833 
58004.0839 


Best OMEGA corresponding to the computed grid-point heights 


z(15) 
0.882 
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is too large by a factor of 3. The number of iterations rises from 
3 at the second grid point to 10 at the sixth grid point, jumps to 60 
at the eighth grid point and would never converge at IKM = 9 because 
there is an oscillation between two stable values at this grid point, 
Examining the values of OMHGA at various stages of iteration, we obsere 
ve that a different value for OMHGA is needed depending on whether we 
start from below the actual value or from above. This phenomenon peaks 
at IKM where the iteration process breaks down if we used OMEGA = 1, 
At IKM it converges to one value above the actual value for z and to 
another value below the actual value for z. This points out that there 
is no guarantee of an unique answer with this iteration process. In 
fact thare are three solutions to the equation for ae one using the 
lowest value for z in the linear term and the highest one in the loga- 
rithmic term, another one using the lowest value for z in the logarith- 
mic term and the highest one in the linear term, and finally one using 
the actual value for z in both terms. Obviously we are interested 
only in the last solution. The only iteration process which would 
give an unique answer all the time is the one onY, but this is 
generally a fairly slow process and difficult to implement. Now we 
will try to derive an analytical function which will give the best 
OMSGA at all levels below IKM. The matching can be done with the 
help of the following rules, 
(a) exact match of OMEGA = 1 at k = 1. 
(b) Match the lowest value for OMEGA at IKM, i.e. OMHGA = 0.2695. 
(ce) Assume an exponential variation of OMBGA in terms of k, 

Therefore, we will use an expression like, 


OMHGA(k) = 1 + af teoxp[(ke1)/(tav1)}"} (A.12) 
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where A and n are constants to be determined by fitting the actual 
values found for OMEGA at various grid points. Good values for these 
constants were found to be, 

A = 0.425 (A. 13) 
n = 3 (A.14) 

It is obvious that this formula will not be perfect for all 
roughness heights because of the arbitrariness in the choice of IKM 
and the discreteness of the levels which affect the efficiency of 
(A.5) at IKM. However, it is clear that an uniform value for OMBGA 
at all levels is not acceptable, A value of 1 would not give convere- 
gence to the desired value at IKM, whereas a value of 0.3 would 
increase considerably the number of iterations required at the grid 
points where convergence is rapid with OMBGA = 1. Formula (A,12) 
seems a fair compromise between simplicity, generality and speed of 


convergence. 


A.4.2 Analysis of the Iterative Process Near and Above IKM. 

It takes many more iterations near IKM than much above 
that level. We also note that the initialization improves as the 
level of the grid point increases. In this case the best OMEGA 
varies more gently above IKM than below it. Therefore, we try to 
fit an analytical expression of the forn, 

OMEGA(k) = 1 = B (NS=k)/(NS-IKP) * (A,15) 
whore, NS = top of the boundary layer = 40 
IKP = IKM + 1 


exponent determined to be about 5 


Ma) 
u 


constant determined to be about 0.32, 
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A.5 Iteration Process on & ° 

In the cases where (A.5) nor (A.6) have not converged thin 
the specified number of iterations, we can complete the process by the 
iteration on & which uses (Aoi) in the forn, 

yt = actin[(2tt2g)/2gl/ko + 2t/2 ot — (Ae16) 

If we define DIFF Ue - Y - Then, when DIFF is positive, the 
implication is that 2* is too high and similarly, when DIFF is negative, 
z* is too low. Therefore, the next guess will be, 

2" = 2" = f(DIFF) (A.18) 
where f(DIFF) is a positive function involving DIFF. This function 
has to be determined empirically and will probably vary from grid 


point to grid point. 
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APPENDIX B 


STEADY-STATE ATMOSPHERIC EQUATIONS 


Bel Equationss First Fc -m, 


The basic equations in log-linear coordinates are given in 


Section 3.3.1 and simplify to, 


lees oe ak Tages vege) ) = £(V-V,) (B.1) 
Ba Nis <5) 4) 
dy? irs Hoge < See = £(u,-U) (B.2) 
~ 1215/7, p(2¥ Se! 472 ane nie ay 
A:45((0  (2%)?)8 oe pas: ee) 26K, \8¥) "(34)" 
- 5) | (B.3) 


The diffusion operator )) ( ) is defined in (3.30) and used 
in the equation for K,; 


ah f (14304 )0340.23x4A4 1. (02)) 
02 +6 a5)? ae (24) 


The other coefficients of diffusivity are not used in the 


(B.4) 


computations but can be evaluated diagnostically, 


Ah (e? + 4x0.23 A, 1 £ (0%)-6Ay Qe LY] (24) ) 
e° (B.5) 
The other turbulent moments involving 9 or q vanish, 
Q2 = 6q =q2 = 0 (B.6) 
The mean mixing ratio and the mean potential temperature are 


assumed constant throughout the boundary layer and their actual values 
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are determined by their values at the surface of the earth. 


Be2 Mixing Length. 
The formula for the mixing length is, 
k (2tz,) 


f= =~ ol 2a9) (B.7) 
ioe 


where 9 . is defined by, 
(- ze dz 
Q. = Luecics (B.8) 


9 2 a2 


These integrals are transformed into finite summations by, 


i: Ze dz =|". @ az dy = Sa z(k) e(k)/DELZ(k)  (B.9) 
° ° dy k=1 


F z dz a Z om dy= "a 2(k)/DELZ(k) (B.10) 


Now we look at the vertical derivative of y « in the surface 
layer | varies exponentially in Y -coordinates. A simple centered 


finite-difference approximation like, 
ie me, ~ J (ie #1) 
ay k 2 
would have a large truncation error. Therefore, we will use a more 


exact formulation obtained by taking the derivative of (B.7) with 


respect to ¥/ : 


QV 
— 


= DELZZ(k)/)*(k) (Be11) 


— 
~< 


where yy }° a” 


DELZZ(k) =< 34]2 2 (Be12) 
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B.3. Vertical Derivative of Kye 
As we know, in the surface layer K, behaves like, 
fees uk (2420) (B.13) 


The grid is predominantly logarithmic in that region, which 


implies, 
Y % In[(2t25)/29] (B.14) 
Therefore, | 
Ko exp(Y) and (B.15) 
7 x exp(¢ ) (Be 16) 


Therefore, a simple finiteedifference approximation to the 


vertical derivative of Kw 


K (kt+i)@K (ke1) 
0) Bete eee ene = (B.17) 


will be greatly in error in the lower boundary layer. The following 


transformation corrects the problem: 


do Km _ a1n Kn 8 
7 Ka oT, (B.18) 
This is better in the surface layer because 
Ink, ~Y and (B.19) 
In K 
ce m & constant (B.20) 
a7 


The finite-difference approximation in the lower boundary 
layer should be 


eS | K(k) 


oY 2 In| (Ky(k+1)/Ky (ke 1)) (B.20) 
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Bo Finite-Difference Equations. 
Be4+.1 Mean Wind U. 
The appropriate equation for U is (B.1) in which the various 


terms are approximated by, 


i | oe K (kt! K,,(k-1) above the constant-flux ees 
Kf) fe in the constant-flux layer 
2 
A(2)= Sin 24. )745, 2 1 = DINK DELZ2(k)/24K,(k) DEL2Z(k)(B.22) 
Y Zz 327 
A(3) = Ky(%)* = K (ie) DELZ2(k) (Be) 
A(4) = ay = [U(k+1)-U(k-1)] /2 (B.24) 
oe: = U(kt1) = 2 U(k) + U(k=1) (B.25) 


We group all the terms different from U(k), 
DU = [ U(k+t)+0(ke1)} A(3)4A(4) A(2)4F (V(k)=V,) — (Be26) 
Finally, 
U(k) = DU/[2 A(3)] (B.27) 


Bo4.2 Equation for the Mean Wind V. 
The terms A(2) and A(3) will be used also in the equation 


for V. The other terms are 


A(5) = oa = (v(k+1)-V(k-1)] /2 (B.28) 
xy = V(kt+1) = 2 V(k) + V(k-1) (B.29) 
he 


All the terms except V(k) are grouped into, 


pv = [v(Ke1) + V(k=1)] a(3) + A(5) A(2) = F[U(c)-UZ] (B.30) 
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Therefore, 


Vix) = dv/2 a(3)) (B.31) 


Be4.3 Equation for o. 
The equation for e is (B.3) in which we have to evaluate the 


vertical derivative of (el): 


pare er (B.32) 


We will use a simple centered finite-difference approxima- 


tion to the vertical derivative of e such that 


y = [e(kt1) - e(k-1)]/2 (B.33) 
The vertical derivative of f is given by (5.11) and will be 
called DL(k). Therefore, (B.32) can be written as 
DEL = e(k) DL(K) +(x) [olicti)=o(ket}]/2  (Be34) 


Io? 
7 


A(1%) = (DEL DELZ2(k) + e(k) 2 (ic) DEL2Z(k)) 0623/2  (B.35) 


The coefficient of 2 becomes 


32 eA 
The coefficient of Sye is 
A(13) = e(k) A (xk) DELZ2(k) 0.23 (B.36) 


2 
= ora Lacs) AC4)4A(5) a(5)] DELZ2(k) 
(B.37) 


We isolate the terms in e@(k+1) and those in e@(k-1) and 
call AA and AB their respective coefficients, 
AA = A(13) + A(14) (B.38) 
AB = A(13) = A(14) (B.39) 


Gathering all the terms except e*(k), we have 
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A(16) = e@(k+1) AA + *(k-1) AB + 1.2 K,(k) A(8)  (B.39) 
Finally we obtain 


e*(k) = A(16)/A(15) (B.40) 


e(k) = sort[e?(k) } (B41) 


Be4.4 Equation for Ke 


The array B contains most of the commonly-used variables 
which do not need to be computed each time that the program is run, 


In this case we have 


B(28) = (1+3C,) (B42) 
B(61) = 3 ke (B.43) 
B(5) a Ay (B44) 


The diffusion operator applied to e“.is called A(9) where 


—a(9)= 4 as) Dye?) = B(3t) Acie) [ad 02(iett) + AB 0%(ket) = 2 A(13) 


e2(i) | (B45) 
A(10) will represent e4(k) so that 
A(10) = 0#(k) e(k) (B.46) 
Tne denominator becomes 
A(12) = e@(k) + 2 B(61)A2(e) a(8) (B47) 
A(i1) = A(10) B(28) + A(9) (B.48) 


Therefore, K Ck) is expressed as 
K(k) = B(5) f (ke) AC 11)/A( 12) (B.49) 


Bo4.5 Iteration scheme. 


The accelerated Gauss-Seidel iteration process is used for 


U, V, K,, and e. The basie scheme has been described in detail in 


(neve): (8) ot Set remenssn 
. aa) aes ere eS 
ms ane rn : 


* bg ah ed ()%aftaee = bate 


Py a 
P * ae. ae Ay '% 


Ly 7 2 os tah 


caldata bonsivetinnamne oda ee # es a4 / 


Yee D vy bea Bs 


iment sama aie nae 


ie 


u 


(i OT OG ee ae 
(Pa.@) aa ig ue eo taie aid 
cr 8) 7 : a 7 7 an re 
wcvnats CO), dealin er ow battgen reas ascot sth 
if? )A Se (june ll " (roi satel 5 . os . 4 
easy. \ 4d ae 


~ 


\ - diay wa 


(Oe. td enya 
ei TI a 


g 


Cees "EREDAR eh eg eS? 
ee et £ ‘ iT ‘ ont at sue Ge Sy ' iy 
P - ; : - @ 4 


“a a 7 
i ol > u 


he 
(OF (fol he: Seema tice 


182 
Section 3.5. Some refinements were needed in order to prevent divere 
gence which can occur very easily in the uppermost levels. Most of | 
the problems occur when K, drops to a very small value which causes 
the wind to become very large. Therefore, a minimum value much larger 
than the molecular diffusivity has to be specified near the top of the 
boundary layer. This can be done by specifying a reasonable value 
for Ka at the top of the boundary layer, say of the order of 100 em/ 
Sec, and devising an analytical function for the minimum value of K, 
at the uppermost grid points. In the lower boundary layer the acy eant 
lar value (about 0.15 em“/sec) ean be used as the minimum value for 
K,- When the top of the boundary layer is fixed at 2 or 3 km it is 
necessary generally to compute U and V in double precision in order 
to get a good value for K,. However, K,, e and f have enough 
accuracy in single precision, 

The next problem is to determine when the iteration procee 
dure has been completed. This is done by comparing two successive 
iterative values of K,. We could use a criteria expressing the mi- 
ninum percentage that the difference between two successive iterative 
values must meet in order to have convergence. In our case this is 
not desirable because we are more likely to obtain much greater accue 
racy in the lower boundary layer than near the top of the boundary 
layer. If we used a percentage we would have to choose between either 
a fairly large percentage with a reasonable number of iterations or 
a small percentage and an unnecessarily large number of iterations. 
Therefore, it is better to use an absolute number of the order of 
1 to 10 em2/sec as a criterion which would insure high accuracy in 


the region where it is desired, that is, in the lower and middle 
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boundary layer. We will probably have less accuracy in K,, near the 
top of the boundary layer but this has only a very small influence on 


the wind profile. 


B.5 Finite-Difference Equation for K, and Kose 


As K, and K, are required in the time-dependent model they 
will be computed diagnostically from (B.5) which is transformed into 
Ke(k) = Ky) ={B(5) R(x) [acto) + a(9}P/o2(k) ——(B,50) 


Be6 Finite-Difference Equations: Second Forme 

The diffusion terms in the equations for U and V can be 
written in two forms. The first form uses (3.26) and (3.27) and 
has been discussed in the previous sections. The second form uses 
(3.26") and (3.27*), and it has the advantage of giving a reasonable 
profile of not only the mean wind but also of the stress, The second 
form should be used if the set (3.26') to (3.29") is used in the 
time-dependent model. The finite-difference expressions are much 
more stable near the top of the boundary layer when we use the second 
form. This permits a greater OMEGA near the top of the boundary layer 


and insures a faster convergence. 


Be6e1 Equation for U. 
The diffusion term in (3.26') can be expanded as, 


Suns > oy 
Lay Say 3D Se lt 3F SP gr 32 30), k ~2] 32 
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the surface layer 


Ka at om ( Ju,)/ (ac Q) = constant (Bo 52) 


Therefore, it seems justified to assume the following arith- 


metic average 
XKMP = (K,, Sty as $[k,(eti) DELZ(cH1 )4K,(k) DELZ()] (Bs 53) 


The vertical derivative of U is approximated by a finite. 


difference equation centered at kt 


ar ieag 7 DCletd) = UC) (Bo 5+) 


The stress is evaluated at ke} by the same procedure, 


XKMM = (K,, St) os aK (ket) DELZ(Ko1)4K,(k) DELZ(k)] (Be 55) 


_U send 7 CK) =U(ket) (B56) 


2 
oY 


The coefficient of U(k) is 


A(3) = XKMP + XKMM (3.57) 
We group all the terns except U(k) in 
DU = XKMP U(k+1) + XKMM U(kei) + F (V(k) - ve (B.58) 


Finally, 
U(k) = DU/A(3) | - (Be59) 


B.6e2 Equation for V. 
The diffusion term in (3.27") is similar to the one in (3,.26') 
with U being replaced by Ve Therefore, it is sufficient to mention 


the resulting terms. All the terms except V(k) are grouped in 


DV = XKMP V(kt+i) + XKMM V(ke1) +F (U, = U(k)) (B60) 
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V(k) = DV/A(3) 


(B.61) 
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APPENDIX C 


TIME-DEPENDENT ATMOSPHERIC EQUATIONS 


C.1 Equations. 
The basic equations for the timeedependent model in &/-C00re 


dinates are given in Section 3.3.1 and will not be repeated here. 

The equation for the mixing length 5 and for its derivative with rese 
pect to Y have been derived in Section B.2. The discussion in Section 
B.3 about the vertical derivative of K, in the steady-state model is 
also valid for the time-dependent model. However, in this case the 
thickness of the surface layer will vary according to the diurnal 
cycle and may be fairly small during the nocturnal inversion. There- 
fore, we should use (B.20) only up to the level expected to represent 


the top of the nearly-constant~-stress layer during the nightime, 


C.2 Notation. 

In our imolicit finite-difference scheme all of the variables 
in the equations which are not differentiated with respect to time 
have to be evaluated at the two time steps t, and t,t At. This would 
imply repeating two nearly identical terms differing only in their time 
4ndex. In our disevssion we will outline how the calculations are 
done and write out these terms with a general time index t. The fole- 


lowing arrays are frequently useds 


i 


A = expressions involving the atmospheric variables, 


B = various constants, 


expressions related to the computation of the surface tem 
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perature, 
D = expressions involving the atmospheric variables at the 
past time step t,. 
The variables starting with the letter D and containing at 
least three letters indicate a finite-difference approximation to a 
derivative, Termination with { implies the first derivative and 


termination with 2 implies the second derivative. For example, 


DUi = first derivative of U 


second derivative of () ° 


DT2 

The name of the variable being differentiated is found 
between D and 1 or 2. As the symbol @is not generally available 
in FORTRAN, T has been substituted for it. 

As pointed out in Section 3.3.2, the diffusicn term ean be 
expressed in two forms in Y-coordinates. Equations (326) to (3.29) 
involve differentiation-of the coefficient of eddy diffusivity and 
the first and second derivatives of the mean quantities, This form 
leads to computational problems which are described in Section C.5, 
Equations (3.26") to (3.29*) involve differentiation of the flux of 
the mean quantity with respect to Ys This second form gives a realise 
tie variation of the flux in the vertical and is computationally 


more stable, Section C.3 uses the first form whereas Section C.4 


uses the second form. 


C.3 Finite-Difference Equations: First Form. 


Ce3ei Equation for U. 
The basic equation is (3.26). The vertical derivative of 


K,, is evaluated as, 
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3(K, (k+1t)=K,( ke1,t)) above the constant-flux layer 


ae Da = Re as (Cot) 
oY : Ka( kt) ln aeeay in the constant-flux layer 
Other terms involved in the diffusion term are, 
a(2)|, = (un (24 14 +, oY 1 
xy | 
at et DELZ2(k) + Putty DEL22(k)) (C.2) 
ACL) |, = (Sh)? = K (kt) DELZ2(k) (c.3) 


The ratio of these two coefficients is, 
AA , = A(2)/a(1)| $ (C.4) 
The geostrophie and Coriolis terms become, 
= £(VeV,) = 2f(Vksty)W(kstotAt)=2V,)  (Ce5) 


The first and second derivative terms are, 


pvt i 2 oni = U(k+1,t)=U(ke1,t) (c.6) 
2 
pu2| = LU = Uke t)-U(keyt)qU(ke1,t)-U(kyt) (Co?) 
. 0 Y% 


The diffusion term at te is evaluated as, 


D(ky1) = (v2, + An | put], (c.8) 


If we exclude U(k,t,+At) in the diffusion term at t,+At, 


we are left with, 


Dlg bee (Ulett +O t)4U(Ketst +b t) HAA Don), soe) (c.9) 
The coefficient of U(k,t,+At) is, 
= 1/Ot + A(1) | (C.10) 


Therefore, the expression for U(k,t,+At) is, 


U(k, t+ At)=4U(k,t,)/ A t+FAC+A(1) [[pv+D(x,1)) /2}} /CST (C.11) 
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C.3e2 Equation for V. 
The equation for V is very similar to the one for U. There. 


fore, the coefficients of the diffusion term involved in (3.26) are 
the same ones as those involved in (3e25)» The first and second 


derivative terms are: 


Dv, = 2 SS = V(k+1,t)-V(ke1,t) (C212) 
2 
Dv2 I = 3 = V(k+1,t)=V(k,t)+V(ke1,t)-V(k,t) (C.13) 


The diffusion term at t, becomes, 


D(k,2) = (ove |, + AA pvt], ) (C.14) 
fe) 0 


bo 


A part of the diffusion term at t,+ At is, 


DV trot Viktt st tA t)+V(ket yt tA t)+(AA bv1) ) (C215) 


tp tAt 


The geostrophic and Coriolis terms are written as, 


FAC = £(U,-U) = 4t{avg-[V(kyt,)*V(kyt,+at)] | (C.16) 


Finally, 


V(k,t,+5 t)={V(Kst,)/2 ttFAC+A( 1) {([pv4D(1,2)) /2)\ [OST (0.17) 


C.303 Equation for®. 
The derivative of K, with respect to Y is, 


3(K, (k+1,t)=K,(ke1,t)) above the constant-flux layer 


dK (C.18) 
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=DK 
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t | = 
MM), K, (kt) - K, (k+1,t) 
2 Ky (ie~1,¢) 


The coefficients involved in the diffusion term at t, are, 


) in the eonstant-flux layer 


K ; 2 
no, = St GO? +n 53) : 


4 (D«TA|, DELZ2(k) + K,(kyt) DEL2Z(k)) (c.19) 
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(2)? = K,(k,t) DELZ2(k) (C.20) 


The first and second derivatives of @ with respect to Y are, 


A(5) , = 


re 2 - ® -@(x+1,t)-Qike1,t) — (ce24) 


2 =a =O(K+t, ae t)+Q(k-1,t) (x,t) (0.22) 


The diffusion term evaluated at te is, 


DT2 


D(k,3) = (D2 


to ACO)/A(5)DT4) | 4) (C.23) 


Part of the diffusion term at t,+At can be written as, 
DOT | +A 2 Olktt, t+ 0t)+Qe-1,t, +2 t)+{aCor/aCsord) |. nt 
(C.24) 
The coefficient of @W(k,t,+ht) 4s, 
CST2 = 1/At + A(5) (C225) 


Thorefore, 


©® (kt +Ot)= {e (kyt,)/A t+A(5) [[DDT+D(k, 3) /2]§ /este (Co 26) 


C.3.4 Equation for Q. 
The basic equation is (3.29). The vertical derivative of 


ey is, 

Sx | 2(Ky(kt1,t)~K, (ke1,t)) above the constant-flux layer 
Zee =DEWIE = (C.27) 
vy lt , Ky(kt) | sake t) 


a ITED. in the constant-flux layer 


Ghee ike 
The peer of the diffusion terms are, 


ROR <4 ee Pe (2k? +x, 2%, = 


oY on 2 it 
4 (Diwi|, DELZ2(K) + Ky(lest) DBL2Z(e)) (c.28) 
A(8) = K, (24)*= K(eyt) DELZ2(K) (c.29) 


The first and second derivative terms become, 
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Sa ee ia 
rail, = 2 oy" Q(kti,t) = Q(k-1,t) (C.30) 


2 
m2 . = Tu HO(K+1,t)-O(kyt)4Q(Ka1,t)-Q(k,t) (C031) 
The diffusion term at t is, 


D(ky4) = (onal, x(a(9)/a¢@y002)|, (c.32) 


The diffusion term at t,t At can be written as follows when 


we exclude Q(k,t +A 3) 


DDQ tot pp eC kttotyt D tka 1, t,t ht)+(A(9)/A(8)DQ1) tot At 
(C.33) 
The coefficient of Q(k,t + ht) is, 
CST3 = 1/at + A(8) (C.34) 


Finally, 
Q(kst + At)=[O(k,t,)/B tans) [(oDa+D(Ke,4)] /2){ /ost3 (c.35) 


G.3-5 Equation for 6. 
Equation (3.31) 4s the equation for e*. We define some of 


the terms which are frequently used in the equations. The various 


derivatives of the mean quantities are called, 


are 
a21)| ,=2 (PE =2((2 2)? 28) = A (pu1)?x0v1)?) |, DatZ2(1) 
(C.36) 
bral, = 2 als) = 2 oti, DELZ(k) (C.37) 
pon}, = 2 22 = 2 pai}, DalZ(k) (0.38) 


The covariances representing the heat flux and the flux of 


water vapor are, respectively, 


2we => K,(kyt) DTA (C239) 


we, 
wal, = 2 Wa =- K,(k,yt) DOA (C.40) 
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These two covariances can be combined through the concept of 


virtual temperature, 


A(12)| ,= 2 ey = 7 Oy = B(21) WE + B(20) Ww (C.41) 
where, B(21) * ¢/T, (C42) 
B(20) = 0.6078 g (C.43) 
We define, 
ap|.= 02, = elkst) A (cyt) (Cals) 
The coefficients of the diffusion terms are, 
A(13)|,= 22:22 of se B(22) ar], DsL22(K) (C445) 
A(14) | -#( 09], 24s es (22)? ) = ne 
${ ar peL2z(k)+/ ) (k,t) Dal , to(k,t) DL(c)} DEL2Z(k) } 
where, va|,.= 32 Ty = 4(e(k+1,t)~e(k-1,t)) (c.47) 
DL(k) = wi = DELZ2Z(k)/ ) Ce, t) (C.48) 
B(22) = Ags28 (C449) 


The terms evaluated at Ee are grouped in, 


ty ¥Oyr 21 btm ES abs > 
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@B(24) o(kyt,)/ A (kyt,)) (C.50) 
The terms which are different from e%(kyt tA t) and which 


are evaluated at % ot At, AYO, 


A(16) = ( le § 3 (0) 2)} 421K BY Tr w0y)) 


trot 


B(22) (A(13) DB2+A(14) DE1) + K,(k,t) AC11) + AC12) (C51) 
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Finally, 


Bat = oP (yt tht) = (D(iy5)+A(16))/ [2(1/ A t#AC13)] 4g 9g )4B(24) 


o(k,t +A t)/A (k,ytot 4t)] (C.52) 


o(kstyt At) = (6%(k,t + at))@=sort(z2t) — (¢.53) 


C.306 Computation of 9°, 


Equation (3.32) is used. Ths coefficients of the diffusion 
term are similar to A(13) and A(14) which were derived for e*. The 
only difference is that the constant B(22) has been replaced by 0.23. 
Therefore, 


0.6 A(13) 


2 (C254) 


mands 
8 


A*(14) 


Nia 


0.6 A(14) (C.55) 


e* 


The first and second derivatives are rospectively, 


we Re — 
pi2i|, = 2 (20) = Q°(ie+1t) aiortuan,t) (C.56) 


cb 


25 — vs — —s 
DT22 SGnoe O° (ett ,t)~07(,t) +0 (ko 1,6) 0-(k,t) (c.57) 


A term which is encountered in 0%, eq and q*, is, 


aot) eects A teofll Ba} = 2/ D t-B(26) (ky t,)/1(k,t,)(C258) 


ts 


where, B(26) = 2/B, (C.59) 
The terms evaluated at the past time step t, are, 

D(k,6) = (A'(13) DT22+A'(14) DI21-WT pan/2)|, + O2( kat.) A(24) |, 

(C.60) 


:?) 


The coefficient of 02(kt,t Tale ep oe Eo 
A(18) = 2(1/A t4A(13))+ B(26) o(kt,+At)/A (este Dt) — (C461) 


The terms different from §2(kyt + At) at t tAt, are, 
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A(19) = [A(13) Dr22 + A(14) Dr2t - wr paa/a| + nt 6er82) 
3° 


Finally, 


Q7(yt,+ Ot) = (AC19) + D(k,6))/A(18) (C463) 


C.3e7 Computation of q?. 
The reference equation is (3.33) whose formulation is 


similar to (3.32). Therefore, the coefficients A*(13), A'(14), A(18) 


and A(21) remain unchanged. The first and second derivative terms are, . 


=i) 12h iad 
m2i | = 2 = = q@(k+1,t) - q2(k-1,t) (c. 6) 


DQ22 fe * saz = q2(k+1,t)=-q2(k,t)4q2(ke1,t)-42(k,t) (C.65) 


The terms which are evaluated at the past time step, are, 
D(ky7) = (A"(13) Doz +.A'(14) DQ2i = WO DQA/2)|, + a*(k,t,) A(21) 
° 
(C.66) 


We group together all the terms evaluated at t+ At, except 


q*(k,t +At). 
A(22) =A°(13) DQ2 + At(414) DQ21 -wO DQA/2 (c.67) 


Therefore, 


a(t) = (A(22) + D(k,7})/A(18) (c.68) 


C.3.8 Equation for 64. 
Equation (3.34) is used for Ode Again, the coefficients 


A'(13), A'(i14), A(18) and A(21) are the same ones as derived for Q%- 


The first and second derivatives are, 
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rat], = 2 Se = @q(kti,t) - oq(ke1,t) (C.69) 
pre? |. Te =6q(k+1,t)- 6 q(kyt)+oq(k-1,t)-04(k,t) 
(C70) 
The terms evaluated at t, are, 
D(k,8) = (A'(13) DIQ2 + A°(14) DIQ1 -(WQ DTA + WE DQA)/4)|, + 
Balkst,) A(21)], (C71) 


0 
The term at t,+At which is equivalent to D(k,8), is, 
A(22) =(A"(13) DDQ2 + A'(14) DTQ1 (WO DEAWT DOA)/)|, ag 
i@) 
(C.72) 


The final result is that, 


Oa(kyt,) = (A(22) + D(k,8))/A( 18) (C.73) 


Co309 Equation for K,. 


Equation (3.35) is the equation for K. We will drop the 
ime index for the next three sections because the coefficients of 
diffusivity Ka» K, and K,, are evaluated diagnostically at the time 
step t,t At. Some terms are common to all three coefficients of eddy 


diffusivity and ars, 


AG = @3 = 0%(k) e(k) (C.74) 
z 
a(23) =4 ay) 2. (of S27) (0.75) 


When we use the results which we have obtained in the pre- 
vious sections, wo can rewrite (C.75) as, 
A(23) = B(31) (A"(13) DB2 = 2 e2(k) + A'(14) DE1) {(K)(C.76) 
where, B(31) = 4 Ay (C.77) 
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Also, 
an = B(5) Q(e) (c.78) 
where, B(5) = A, (C.79) 
The numerator in (3.35) is expanded as, 


A(24) = AR (B(28) AG + A(23) + (B(43) WE + B(45) wo) V6) (c.80) 


where, B(28) = 1 - 3Cy (c.81) 
21 Al € 8 

B( 43) 2 FE (C.82) 

B(45) = Aish 0.6078 g (C563) 


The denominator is, 


A(25) = o°(kt) + B(47) (AC11) + 1.5 (B(21) DTA + B(20) DaA)) } 2(K) 


(Cc. 84) 

wheres B47?) = 3 AS (C.85) 
Finally, | 

Kik) = A¢24)/A(25) (C.86) 


Equation (3.36) is appropriate for K,- A term common to 


both K, and K, is, 


2 male 
AX = od +4 A, 0.23 | 2- (ol g.2 )-64,) kK, Fes 
= AG + A(23) = B(36) A(t) Ky(ie) AC14) (c.87) 
where, / B(50) 4.5 Ay (C88) 


The numerator in (3.36) is, = 
A(26) AR = AR (AX+43(46) Wa 4 (ie)m6 4 (ic)(B(24) O7(1e)4B(20) Ga(Ie))/DTA) 
(c.89) 


where, B(46) = 6 A, 0.6078 g (C.90) 
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where, 


C.3-i1 


The denominator becomes, 
A(27) = e2(k) + B(58) DIA L2(k) 
B(58) = 60% @/T, 
Therefore, 


K,(k) = AR A(26)/A(27) 


Equation for Kise 


Equation (3.37) is used for K. The numerator can be 


expressed as, 


2(28) an = ar{ Ax+B(63) { (ie) wr-6 0(K) [ B(21) Ba(k)+B(20) 42(x)) /oaa] 


where, 


B(63) = 6A, g/T, 
The denominator is, 
A(29) = e2(k) + B(33) 17(k) Dan 
Finally, 
K(k) = AR A(28)/A(29) 


(C.91) 
(C.92) 


(C.93) 


(C.94) 
(C.95) 


(C96) 
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C4 Finite-Difference Equations: Second Form. 


Co4.1 Equation for U, 
The diffusion term in (3.26") is expanded as, 


2) 
eee 
(C97) 


In order to evaluate the stress at k + $, we recall that in 


P) aU Q 
$5 Spree rs OSH Srey re er 


the surface layer, 
Ky 22 x (2u,)/(ACL) = constant (c.98) 


Therefore, it seems justified to assume the following arithe 
metic average, 
2 XKMP| = (Kn eee gwe(Ky(kt1,t) DELZ(k+1) + KA(k,t) DELZ(I)) 
9 
(C.99) 
The vertical derivative of U is approximated by a finite- 


difference equation centered at kt}, 


dU a f 
y earns U(k+1,t) = U(k,t) (C100) 


The stress is evaluated at k-$ by the same procedure, 


seedy? [Ka (Ketot) DBLZ(lent)4Kq (cot) DBLZ(K)] 
: (C.101) 


2 xa, = (KX, <2.) 


aU = = 
a Tleatst UCk,t) = UCkK-14t) (C5102) 


We group all the results and obtain, 


DELZ (Ke) | U(ictd yt) XKMP-U(k,t) (XKMMAKKMP )+U(k~1,t) 


” 9U 

a (Ky 5) af 
XKMM | (es103) 

We will use the same finite-cifference scheme as described 


in Section 3.4.2. As a result, we group the expressions evaluated at 
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the past time step in D(k,1), 


D(ky1) = DELZ(kX[U(k+H,t,)-U(k t,)] XP], #[U(k-tyty)-U(k to)] - 


xan], } + FLV(kyt9)/2 -Vy]_ (C.104) 
The coefficient of U(k,t,+At) is, 


cst = 1/ A tH(XKHPHKKM) | 5 4 yt DELZ(k) (C2105) 


Finally, 
Ulkstot At) “{[XOP|, ane Ukr t+ At) sme 4 p¢ Ulket,t.+at)] 
DELZ(k)4D(ky1)-F V(k,t,+ At)/2 + Ukyty)/At} /cst (C4106) 
The mathematical expressions could have been condensed but 
' the condensed expressions would not minimize the round-off errors as 
well as the expanded version of the equations. In Section C.6 we will 
discuss in some detail the problem of round-off errors in the determi- 


nation of the profile of the mean quantities. 


Co4¥.2 Equation for V. 
The diffusion term in (3.27') is similar to the one in (3.26') 


with U being replaced by Ve. Therefore, it is sufficient to mention 
the resulting terms. The expressions evaluated at t, are, 
D(k,2) = DELZ(K){ [V(kt1yto)-V( k oty)] XKMP |, +[V(k=1,ty)=V( koto) 


xrcoat| | + F[U, = U(k,t,)/2 | (c.107) 


Consequently, 
G 


V(kstotA t) ={[Hor |, +4 t V(t totAt) + KKM ¢ + qt V(k-1 stot t)] 


DELZ(k)4D(ks2)-F Ulkst y+ At)/2 + Vikst,)/ at{/cst (C.108) 
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C.4.3 Equation for(H). 
The diffusion term involved in (3.28") is similar to the one 
in (3.26") if we replace U by i), and K, by Kye Therefore, parts of 


the heat-flux terms are, 


2 xem |, = (Ke SO] gg = Ele (eetyt) DEEZ (et) + Ky(Keyt) DELZ(K)] 
(C.109) 


oe 


2 XKTP| .=(K, a De 3[K, (kt+1,t) DELZ(k+1) + K,(k,t) DELZ(k)| 


k+3, 
(C.110) 


When we use the centered finite-difference approximation 
to the e/-derivative of ), we obtain, 
2. x ) ke = Dai2(k)[@r+tst) xx0P]_ ~ Olxst)(xxrPKr)| 4 
+ @(k-1,t) xKrm|y J (C.411) 
The expressions computed at the past time step are. grouped 
into, 


D(k,3) = pauz(y) § [Occtt,t,)- Dirt,)] x0 | 4{Q(k-1,t,) - 
@ (sto) ] xeTH| ¢,, } (C.112) 


The coefficient of G)(k,t,+At) is, 
csT2 = 1/ qt HXKTOXETP)| to at DELZ(k) (c.113) 
Finally, 
@ (kytotAt) “Osta, {HP he ent (A(ct1,to+0t)/G(kst,)) + 


XETM|g 404 (leet tot t)/ Olkst,))] DELZ(k)}/cst2} — (¢61444) 


Co4.4 Equation for Q. 


Equation (3.29') is similar to (3.28') if we replace K, by 


Kuro and (f) by Q. Therefore, parts of the water-vapor-flux terms are, 
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2 XKWP|, = (K, ore = $(K,(e+1,t) DELZ(k+1) + K,(ieyt) DELZ(k)] 
Cis) 


ag 7 SK (eet ,t) DELZ(ka1) + K,(kyt) DELZ(K) ] 
(C.116) 


The diffusion term can be expressed as, 
BY 3 sae eS 
$ (25K, 28d], = patz(i)[acer1,t) XeWP] .~o(%e,t) (KKWPAXEWM) |, 
+ Q(ke1,t) XKWM|, | (Co117) 


in which the centered finite-difference approximation has been used to 
evaluate the Ya-derivative of Q. The expressions computed at t, are, 


D(k,4) = DBLZ (ke) [AC e+1 sy )=0C es t)] KKH] +[Q(ke1,t, )=O(k,t,)] 


XKWM| { (C.118) 
° 
The coefficient of Q(ty+At) is, 
CST3 = 1/At+ (XKWP-+KKW) | DELZ(k) (c.149) 
re) 
Finally, 


Q(k,t,+ At) ={ [xxwr eta QCA ty AE) RIG] _ Uke1,t,+at)] 


DELZ(k) + D(k,4) + Q(kyt.)/d t }/CST3 (C.120) 


Co%.5 Other Equations. 


The other equations are left unchanged and can be found 


from Section C.3.5 to Section C.3.11, 
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C.5 Computational Problems. 


In this section we will deal with the computational problems 
associated with the first form of the diffusion term. We will discuss 
only the problems encountered with (3.28) for the following reasons: 
1) The computation of the temperature profile is the most sensitive 
to these problems. 

2) The equation for @) is Simpler to analyse than the ones for U and 
V because it does not involve Coriolis or pressure~gradient forces. 

3) The derived conelusions from the analysis of Hare applicable to 
the other 3 main quantities because of the similitude of the equations, 

Therefore, it 1s sufficient to analyse the behavior of 
the equation for () » Another simplifying assumption in the analysis 
is the use of the forward-in-time finite-difference approximation 
to (3.28) instead of the finite-difference approximation used in our 


equations and described in section 3.4.2. Therefore, - 


oA. 
Nes f (C.121) 


will be approximated in our discussion by, 

A(t,+ At) = A(t.) +At £(t,+At) (C.122) 

The reasons for that simplification are: 

1) The finite-difference equation becomes simpler as it involves 
roughly only half the number of terms required by the other time. 
finite-difference approximation. 
2) The conclusions which will be drawn from our discussions are almost 
immediately applicable to the actual finite-difference equations used 


in the program if we look into the effect caused by two consecutive 
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tims steps. 


Ce5el. Formation of Artificial Inversionss Case of no Real Inversion. 
The inversion can be generated only by the earth's surface 
in our model because the temperature at the top of the boundary layer 
is kept fixed at all times. The top of the inversion will be carried 
upward with the help of turbulent diffusion. The diurnal temperature 
cycle has only one maximum and one minimum per day because we neglect 
advection and cloud formation. The model is started from neutral 
stability at the time of maximum temperature. Therefore, the earth's 
surface will cool down for the first few hours creating a stable 
layer throughout the atmosphere without any inversion. If the model 
develops any inversion, this will be caused by numerical problems. 
The model ran smoothly until we obtained a region where Ky was varying 
rapidly due to the increasing stability, and this caused the formation 
of an unstable layer on top of the surface-based stable layer. This 
behavior caused a deterioration of the model and ultimately caused 
divergence of the results. The finite-difference version of (3.28) 


‘becomes, when we neglect the time-derivative term, 


Ox) Z ONGD Mee ACO Gert) + (AC) = A(2))@) (cx) (C123) 


K,(k) DELZ2(k), associated with the second derivative, 


where, A(1) 


2 K ie 
and A(2) = K, i tirerer iste » associated with the first 
g, rs) 


derivative of @with respect to Y. 
If we restrict K, to non-negative values, then A(1) 
4s always positive. However,A(2) can be of either sign depending on 


the sign of the Y~derivative of Ky. Some possibilities are, 
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Case 13 A(2) = A(1) 


In this case we haves ()(k) = (ict) 


Case 2s A(2) =0 

Now, we obtain: G(x) = Overt) +9 Cent) 
CED A(2) = ~ a(1) 

This implies: &) («) = (4) (ke1) 


Therefore, for 
A(2) € [-A(41),A(1)] 

wa have, cx) € [O(x-1), Oc+1)] 
where € signifies that the term on the left-hand side is within the 
range defined by the terms inside the brackets on the right-hand side. 
Conversely if A(2) is not in the above-defined range, then () (x) is 
computed outside the range [)(k-1),@(k+1)] « This means that the 
diffusion term which normally tends to smooth out any temporature 
gradient, would inerease the temperature gradient if certain conditions 
are met. We will derive the conditions under which such nonphysical 
results are obtained from the finite-difference schemas Mathemati- 


cally, A(2)> | A(1)limplies that, 


2 
3k 2 
be 


Equation (C.124) can be put into a finite-difference formula 
and, for the grid usually used in our simulations, requires a variation 
in Ky of an order of magnitude between two successive grid points. 
Exeessive cooling or warming will be computed whenever condition 
(C,124) is mate This problem can be partly cured by imposing that 


A(2) < \a(1)\ at all time steps. But this is an arbitrary assumption 
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which interferes with the computational scheme, The inclusion of the 
time derivative term adds very little, and the condition A(2) < \A(1)| 
insures now that, 


Dlkst tat) € [lst Oro1,t,+ At), O(k+1,t,+0t)] (c.125) 


C.5e2-e Blocking of the Height of the Daytime Inversion. 


Another problem related to the use of (3.28) appears a few 
hours after sunrise when the inversion height, after rising relati- 
vely well, suddently refuses to move higher than abcut 50m. This 
produces a very nneteels surface layer capped by a very stable layer. 
The causes for that behavior are a complex interaction between Ky» o2 
and @. 

1. Ky depends strongly on the sign of the temperature gradient. A 
small value of K} is obtained in stable cases and a large one in un- 
Stable cases, The difference in sign of the temperature gradient 
means usually a few orders of magnitude difference in the computed 
value of Kyo 

26 G2 should nearly vanish whenever the temperature gradient decreases 
to a very small value. Below the inversion, in the unstable region, 
2 is relatively large. In the stable region above the inversion, 

92 is again important where the temperature gradient is large. There- 
fore, we expect a minimum in 92 at the inversion. However, the value 
of this minimum becomes dominated by the diffusion term and the time- 
dependent terms in (3.32) and no lenger by the temperature gradient 


term, That means that 6% does not vanish in the limit of a very small 


temperature gradient. 
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3. The equation for(H) 4s not solved in terms of the vertical variae 
tion of the heat flux. Consequently, the computed temperature profile 
cannot force alone the lifting of the inversion height, but permits 
very often the existence of a very unstable layer capped off by a very 


stable layer, 


In the equation for Ky the product 0 20 appears, Near the 


inversion, 92 is larger than expected from the a A® 70 and 2@ 


is positive. This gives a very small value to Kes which in turns 
prevents the lifting of the inversion height. In order to lift the 
inversion height we would need/\@ “0 , which causes large values for 
p2 and K,. This would maintain the inversion height at a relatively 
high position. The computational scheme seems to be able to lift 
the inversion height under some conditions which are not always met. 
One way of euring the problem is to neglect the time-dependent and 
the diffusion terms in (3.32) near the inversion height. This allows 


a smaller value for 9“ 3 near the inversion height, which permits 


a larger value for Ky and which ultimately will lead to the lifting of 
the inversion height. However, this method of solving the problem 

can initiate some oscillations caused by the variations in 2 because 
we are using two different formulas to compute 02, Therefore, convere 


gence is not readily insured at each time step. 
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C.6 Refinements to the Finite-Difference Equations. 


Ce6.1 Double Precision, 

Near the tcp of the boundary layer, the mean atmospheric 
variables vary slowly both in the vertical and in time. However, the 
coefficients of eddy diffusivity and the other turbulent quantities 
vary appreciably in that region and their exact value depend to a 
great extent on the gradient of the mean atmospheric variables. There- 
fore, in order to get a good estimate of the turbulent quantities we 
need a fairly accurate evaluation of the mean variables, This implies 
that double precision is generally needed for U, V, Gana Q whereas 
Single precision is probably always sufficient for all of the other 
atmospheric variables. The great inconvenience of double precision is 
that it requires 2 to 4 times the computing time of single precision. 
Consider a specific example in which we assign 7 decimal digits to 
single precision and 14 decimal digits to double precision. We define 
an elementary operation as the operation between two digits. We will 
compare the number of elementary operations required for addition and 
multiplication in single precision with the number of similar opera- 
tions in double precision, We consider two numbers A and B which 
consist of n digits each, 

A = X4 Xp e098 X, : (C5126) 
B= ¥, 99 °t* Ip (4127) 

The addition A+B involves n elementary addition of the type 
X,tyy- If x47, 7 10 we add i to x,_,. There is on the average a 50% 
probability that this additional operation has to be done. Therefore, 


an addition requires 3 n/2 elementary operations. 
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Table 10. Average number of elementary operations for addition and 


multiplication in single and double precision when we use the decimal 


precision | single 
operation 7 digits 
0 


addition af 
multiplication 126 


Table 11. Average number of elementary operations for addition and 


system. 


double 
14% digits 


multiplication in single and double precision when we use the binary 


precision | single double 
operation 7 digits | 14 digits 
addition 9 18 
multiplication 102 399 


system. 
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By similar reasoning we can estimate that an average multi- 
plication will require about [ (n=1)® + n2}3/2 elementary operations. 
Table 10 compares the number of elementary operations required in single 
precison and in double precision for the operations of addition and 
the operations of multiplication in the decimal system. Table 1i does 
the same comparison in the binary system. The most important onelusion 
which can be drawn is that addition requires only about twice as much 
computing im double precison with respect to single precision, whereas 
- multiplication in double precision needs at least 4 times the compu. 
ting time of a single pecision multiplication. iiherecare if we want 
to save computing time, we should try to reduce the number of multie 
plications which are done in double precision. For example the tempe-~ 
rature difference @ (k+1) - ) (ke1) is computed in double precision 
and the resulting value is put in a single precision variable DTi. 
This temperature difference is needed, for example, in the computation 
of some of the turbulent variables. By using DTi we do not lose 
any significant information and we gain by reducing the number of 


computations performed in double precision. 


C.6.2 Minimization of Round-Off Errors. 

The round-off errors ean be very important near the top of 
the boundary layer. During seme numerical experiments these have been 
observed to add up consistently, This effect was most visible on 
the temperature profile which showed an increasingly large tempera- 
ture gradient between the top two grid points. When the computations 
were done in single precision this was enough to prevent the convere 


gence of the model. In double precision the problem is somewhat 
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reduced but problems can arise if the model is run long enough. 
Therefore, it seems important to devise a finiteedifference scheme 
which reduces the possible round-off errors, This is achieved mostly 
by being careful in choosing the order of the operations inside an 
expression. We will take the example of a constant isentropic atmose- 
phere, for which a round-off error will mean a deviation from the conse 
tant value. Ordinarily the computer transforms the decimal numbers 
into hexadecimal numbers. There are a few operations in real number 
arithmetics which should not give any roundeoff errors. These are 
1- The hexadecimal representation of a small real decimal number 
without any decimal part is exact. For example 3. would have an exact 
hexadecimal representation whereas 0.3x102 would not. 

2- The difference A =~ A = 0. exactly. 
3e The division A/A = 1. exactly. 
4. The sum A + 0. =A exactly. 
5- The operations Axi =A and A+1=A exactly. 
6= The addition or subtraction of two small real numbers without any 
decimal part is exact. 

There are other operations which may give round-off errors? 
“4-2 The multiplication between two real numbers A and B, followed by 
the division by B, gives A with a round-off error. (A x B)/B =A + round- 
off error. This round-off error could be prevented by changing the 
order of the operations: A x (B/B) = A x (1) = A without round-off 
error. 
2. The operation (A + A)/A gives 2 with a round-off error. This can 


be cured by (A/A + A/A) = 1. +1. = 2. exactly. 
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We go back to (C.26) in which D(k,3) is identically zero 
for a constant-temperature profile because DTt | and DT2 to are eXaCe 
tly zero according to operation 2 which is the difference between two 
equal numbers. We note that Drs | tot Ot is also identically zero for 
the same reason. We are left with 
@(k,t.+ at) AOx,t,)/ Att A(5) 13 [@ (ttt, + ht)+Q(k-1,t.+ aryft / 

L1/ 0 t+A(5)) 

The operation cet t+ dt)+O(ke1,t,+0+)} may or may not 
give a round-off error in the case of an isentropic atmosphere. For 
example, we suppose a foupedieve accuracy and @= 50103. The real sum 
is 1002.6 which is transformed in four-digit accuracy into 1003. with 
a round-off error of 0.4%. The subsequent division by 2 yields 501.5 
with a round-off error of 0.2. However if = 301.3, the sum is 602.6 
and the divison by 2 gives back 301-3 with no round-off errors. There- 
fore, wo suppose temporarily that there is no round-off errors for that 
operation. We have the ratio 

© £ (W/ At Ay Cy 
1/ bt + AC5 

The operation ® / Dt and A(5) (i) are very likely to have 
round-off errors as well as the subsequent addition. Therefore, that 
formulation is very likely to give rise to round-off errors. However, 


if we rearrange the terms as, 


@ = @ HAE | 


The factor 1/ At + A(5) is common to both the denominator 
and the numerator and will give some answer B with round-off errors 
with respect to the true solution. This is not important as B/B = 4 


without any round-off srror, and @is computed witnout any round~off 
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error. Therefore, the best finite-difference equation isolates 


and factor A(5)- Equation (C.26) is rewritten as 
Mekytor ht) =Ocsty){ (1/a t4a(5) )ffovr/ Olest9)+0(K,3)/ Hest] /2} 
/(1/ot + a(5))} 

We have noted that in actual numerical tests DDT/(H)(k,t,) 
gave exactly 2. with no round-off errors. This implies that there is 
no need to extend the division inside DDT, The rearrangement of the 
terms is effected for all the other finite-difference equations for 


the mean quantities. 


C.6.3 Higher-Order Approximation to the Derivatives. 


In Section 3.4 wa have devised finite-difference approxima= 
tions to the first and second derivatives which were more exact by an 
order of approximation in the Taylor Series. As pointed out 
previously there is no need for the improved and time consuming 
second-order approximation to the derivatives in the regions where 
the atmospheric variables are varying gently, namely in the surface 
layer and in the upper boundary layer. The higher-order approximation 
is optional in the time-dependent model and is used in the region 
determined by the range [m(4) ,m(3)], where mis an integer variable 
used in the program. 

C.6.4 Restriction on the Temperature Profile. 

As noted in Section C.5, (3e28) has a tendency to create 
many stable-unstable transitions. We can devise limits to the tempe- 
rature at each grid point based on the principle that there is no 


internal energy source term and that warming and cooling can only 


» eit: Yo. deimeageetiaet, eft = mt eb he 


we ene diva wt st 


ish [Cd i y ym 
{aherawn’ faa 


ut 4@)\r00 aioe? Lapbaeanee' faudpiy nt te? 


chats Judd setiautipaer, ean 


neh -eotdastpn iaUmPRa HS . 


wt ROA wcciathieldeas Ghai 
oe Xi lea mops. gnu eet seine are Aaa 
St tobi a seth 20A8 path 

new dee a tare oem ods a0 0 

_ ede ea a a nk Te, OAS OF hector setnoxbooaes 
calidon si alk Sees. \yaltony, piieion ova alsa won, 


aheets of vormion! aed ‘asi oil nohecee ode ronan ) 4 >. . 
ieewd eat of sSinth wabreh ra ee. sti ton alg int Rita 


r a wn ut evens fads ai ary bared fale ning isu grind 
“a7 vino, ney pellcos ies. : As mane! emetic sie tact f 
" g -". - Pe FU Ss us “ye 
Wea of ot 6g fie eee 


21g 
come from the other levels. In the program we have called 
D(kys10) =Xk+1,%, =H) (koto) 
which represents the stability at the past time step above the grid 
point Koe At the actual time step, the stability around k, is defined 
by 
SIGN = Mix +1,t,4d t)-Qky-tstotA t) 
From the knowledge of these two variables wae want to impose 
a restriction on @(k,st,+ At). As an example, we suppose that D(k,, 
10) is positive, ener implies stability above the grid point k, at 
the past time step. If SIGN is also positive, that implies that the 
layer around Ky at the actual time step is also stable. The restriction 
to @(k,,t + Ot) is to impose the condition that it must be smaller 
than (D(kot1,t,+At). Only when H\(k,-1,t,) is greater than ()(kgst,) 
is it possible to expect that the restriction interferes with a 
real change in stability, But this case is certainly rare as it would 
require very rapid changes in the surface temperature. We tabulate 
all the possible combinations of D(k,,10) with SIGN and the restrice 
tion imposed on G) (k53t,+ Oe) in Table 12. This restriction on the 
temperature is optional in the time-dependent model. It can be imposed 
for any layer that we wish, it can be done during the computations 
at each grid point or after a complete iteration has been completed, 
Whenever the restriction is imposed, we reset (kot tht) wtniy the 


range [@(kyttot + Ot), @ (kj-tetoth t)] - 
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Table 12. Restriction on((k, st, +At). We defines 


D(k910) = GQ (k,+15t,) -Hirgst,) 
and SIGN = @(ko+1,t +A t)-G) (ko-tst,+ At) 


£ Olkottt,+at) | > Qikg-t,t,+ at) 


{ (x -1,t,+4t) 


D( ky 5 10) 


positive 


negative 
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APPENDIX D 


INFRARED FLUX DUE TO WATER VAPOR 


Del Hypotheses. 
The following hypotheses are made in order to simplify the 


computational task; 

Hypothesis 1: the pressure at the roughness height is constant at all 
times. This implies that the pressure at the top of the boundary layer » 
will vary. 

Hypothesis 2: we neglect the infrared flux emitted by the layer bet- 
ween Zz, and the ground. This layer is fairly thin and the infrared 
flux received at the earth's surface is not modified sensibly by ne. 
glecting the contribution from that layer. 

Hypothesis 3: the eee above the boundary layer is not modelled 

and we neglect the infrared flux coming from that layer. That assump- 
tion is justified by the relatively rapid decrease in the mixing ratio 
with height and by the fact that most of the infrared flux emitted by 
water vapor and received at the ground is emitted in the lowest 100 
meters of the boundary layer. 

Hypothesis 4: only the mean infrared flux will be considered and con- 
sequently we will use only the mean atmospheric variables in our com- 
putations. 

Hypothesis 5: the hydrostatic equation will be used to compute the 
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De2 Computation of Temperature. 
In the computation of the infrared radiation we need the 


temperature which ean be obtained from the definition of potential 


temperature. In index form we have 


K 
nx) = oo [ 2] (D.1) 


where K= R/C, and P, = reference pressure = 1000 mbe 


D.e3 Computation of Pressure. 


We use the hydrostatic approximation 


y) P = @ 
Ste Pg (D.2) 
The perfect gas law is 
P= p R cS (De3) 


The virtual temperature T,, is related to the temperature T 
and to the mixing ratio Q by the approximate expression 


Ty = T (1 + 0.6078 Q) (De) 


Combining these equations with (D1) we obtain 


Kot oP g Po 
Pe De 7 7 HUD H0.6078O) bar 
1 Pee a (D.6) 
K yz RG (170.6078 Q) : 


Equation (D.6) has isolated the pressure and will be used in 


a finite-differenee form 


= pide Po tant) 
P(k+1) = | Cy PU) [O(c+t)+O(k)] {140.3039 19(k+1) 7k J} 


(D.7) 
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De4% Equation for the Infrared Flux. 
The basic equations for the infrared flux are (5.9) and 


(5.10) which we recall here, 


n 
) 
Fy(0) = 5 (oT), ($5) Ou (D.8) 
Nu oa P Q vd (D.9) 
where F (0) = infrared flux received at z, due to water vapor 


n = index of the top level considered 
€ = emissivity of water vapor 

P = average pressure of the layer [mb ] 

Q = average mixing ratio in the layer (g/g) 

g = gravity = 980 em / sec” 

t= average temperature of the layer (°K) 

Au = path length in a given layer (cm) 

O = StefaneBoltzmann constant. 

We will consider the layers between grid points. An average 
quantity in the layer will be taken as the arithmetic mean of the value 
of that quantity at the upper grid point and of its value at the lower 
grid point. For example, the average pressure is, 

P =[P(k+1) + P(k)/ /2 
Therefore, (D.9) is written in finite-difference form as, 
Au =[PCctt)#P(ic)] [Pct )=P(k)) (Q(k+1)2Q(k)] /(4 g) — (D.10) 

The total path length wis defined as tho sum of all the 

path lengths contained in the layers between Ze and the top of the 


layer under consideration, 
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(top) 
we Aux (D.11) 
(z,) 
The slope of the emissivity curve can be derived from 
Table 2, 
my Ory 7 
28 = IV for O<w<i0™ (D. 12) 
: a aa 
26 = 240398 for 10° «ew < 107? (D.13) 
= = 060565 for 10° %w < 1072 (D. 1) 
> = 206 ~ = 
2<. = Oaoet for 10 oe <107! (D.15) 
QE. = 00778 for 10° Ley (D.16) 


Y 
r=] 
x 


The finite-difference equation corresponding to (D.8) is, 


aS x (att € 


D.5 Computation of the Infrared Flux. 

All the quantities can now be computed and the successive 
steps in the computation of the infrared flux are, 
D.5-1 For the Lower Boundarys 
Step 1: Compute from (D.1) the temperature at 256 
D.5.2 For the Other Levels. 
Step 2: the pressure at the grid point k+1 is computed from (s7)s 
Step 3: the temperature at the grid point kt1 is chiubatee from (D.5). 
Step 4: Au is computed from (D.10). 
Step 5: the total path length w is computed from CDaliv. 


Stevo 6: the contribution to the infrared radiation from each layer is 
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computed from (D.12) 


Dee. Integ rated Result. 


Step 7: finally the contribution from each layer to the infrared 


flux is added up to give the total infrared flux received at the 


ground. 
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APPENDIX & 


SOLUTION OF THa SURFACE Ta&MPERATURE EQUATION 


We will use the results obtained in Chapter V and trans- 


form the equations into a finiteedifference form. Equation (5.20) 


becomes, 
EB = KW(2)x(Q(3)-Q(1) )*DELZ(2)/2 (E.1) 
where KW(k) = coefficient of eddy diffusivity for water vapor at 
grid point k Se 1! 
DELZ(k) ® the derivative of Y with respect to z at grid- 
point k 
Q(k) = mixing ratio at grid-point k 
Wquation (5.19) is transformed into 
S = KS(2) (TS(3)-TS(1)) DD2(2)/2 (E.2) 
where KS(1) = molecular coefficient of diffusivity in the soil 


at grid-point 1 = § a pa 

TS(1) = soil temperature at grid-point 1 

DDZ(1) = first derivative of § with respect to 2, at grid- 

point 1. 

The equation for the atmosphericesensible-heat flux involves 
the potential temperature whereas all the other terms in (501) contain 
the temperature. This is why we trensuseraGeea) into a form compati- 


ble with the other terms. The equation for the potential temperature 


is 
1000 K 
@) = 7p Seager (E.3) 
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where k= R/C, = 2/7 


Taking the derivative of the logarithm of (E.3), we have 


ol . olnT 1 
2in@ . dint, Qinp (Bo4) 


OZ 


We use the hydrostatic approximation in (E.4) so that 


9@ _W, ar 
oS) (.5) 


Equation (E.5) becomes, in a finite-difference forn, 
A@= (AT +hz g/C, O/T (E.6) 


The integral form of (5.13) is 


H 
A® = jE, Pry In ( 
3 Ko Uy 


ZtZoe 
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) (B27) 
6 


We combine the last two equations and obtain 


lag LO2)- (4) 
fam Pry fs (ztz5)/z ( M2)= 14) ) awa} & 2(2) (8.8) 


In the last expression we have evaluated the finito-diffe- 
renee equation betwee Zo and the first grid point. We combine the 
outgoing terrestrial infrared radiation with the sky infrared radla- 
tion due to CO> 

RT = E070) (£09) 
where E=(1 =0.18)€, = effective emissivity of the ground. 

The infrared radiation from water vapor (Re) has been come 
puted in Appendix D, and the incoming solar radiation (i) has been 
defined in Chapter V. Therefore, the following equation results if 


we group all the results, 


T(0) + c(2) T(0) + C(41) = 0 (E610) 
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where, Ps cs Ks DDZ(2)/2 + BR 
(2) =—————— 


Eo 


m= pc ky Ux (Ge +@an 
P Pp Po (in (242,)/2,) The) 71 


(IetPt L K,(2) Q(3)-Q(1) daLe(2)/2epB 1(2)422(2) 


hg ae ea: 


Equation (E.10) ean be solved exactly. We will follow the 
steps outlined in the i4th Edition of the CRC Abridged Mathematical 
Tables (1967), pp» 291-293. This requires that we solve first a re- 
solvent cubic equation, whose solution gives the desired solution to 
the quartic equation after a few manipulations. The procedure will be 
simplified by making use of the fact that C(1i) is always negative and 


that C(2) is always a real number. 


E.1i Solution of the Resolvent Cubic Equation. 
The quartic equation to be solved is of the type 


stead to xe tex+d =0 (Bei) 
where a2=b=0 
e = C(2) 70 
d= ECG 10 


The resolvent cubie equation is 
pteytf=0 (B12) 
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1 
Beet. £2. et 3 
p= (-£4(£ +-)4) (B14) 
The three solutions are 
yi - A+B (B15) 
Yo =~ $[(a+B)=(a~B)(-3)4] (3416) 
Yq = = $[(AxB)~(A~B)(-3)? ] (2.17) 


If the disciminant ae +3 ) is 
> 0 wo have one real root and two imaginary roots, 
=) we have three real roots of which at least two are equal, 
£0 we have three unequal real roots. 

In our case 


e > 0 because e is always positive, and 


£e > O because f is real. 
Therefore, the discriminant is always positive and the 
real root is Vye We define the following variables: 

C(3) == 4 £ = C(2) c(2)/2 (8.18) 
c(4) = £2/4 = ¢(3) c(3) (E619) 
C(5) = 0/27 = -[o(1)/4 64/27 (B.20) 

The discriminant is 
DIS = C(4) + C(5) (£21) 


The constants A and B are, therefore, 
o(7) =a =(c(3) + sorr(pts)) 73 — (w.22) 


c(8) = B =[c(3) = sarr(pts)}/3 (#23) 
where SQRT = square root of the function inside the parenthesis, 
Therefore, the solution of the resolvent cubic equation is 


Yi = c(7) + c(9) (E.24) 
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We note that B is always negative but smaller in absolute 


value than A, so that Yi is always positive. 


Ee2 Solution of the Quartic Equation. 


We defins R as 


R= ( a2/4 = b + 1) = SQRT(Y1) (£625) 
Therefore, R will always be real and positive. In this case 


we define D and 5 as 


D = (-R2 = 2 ofR)? (8.26) 


B = (-R2 +2 /R)? (B27) 
We must eliminate all of the non-physical roots. ec = C(2) is 
always positive which implies that D is an imaginary number and is 
rejected for that reason. E is the only real root for which the folloe 


wing inequality must be satisfieds 


-R* +2 ¢/R 7 0 (8.28) 


The two possibie solutions are 


+ a 
x 3 GeiR a auks) (B.29) 
As R and EB are both positive, the only positive real solu- 
tion of the quartic equation is 
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APPENDIX F 


LIST OF SYMBOLS 


kinematic coefficient of viscosity for heat 


also roughness angles angle between the lower boundary 
layer wind and the geostrophic wind 


ratio of the heat capacity at constant pressure over the 
heat capacity at constant volume for air 


Kronecker's delta 

also solar declination 

finite difference operator 

partial derivative operator 

alternating tensor 

also ratio of the gas constant for dry air over the one for 
water vapor 

also effective emissivity of the earth's surface 
emissivity of the earth's surface 

latitude 


ratio of the gas constant over the heat capacity at conse 
tant pressure for the air 


various length scales 

first and second coefficients of viscosity 
molecular coefficient of diffusion for heat 
kinematic coefficient of viscosity 


density of the air and density of the soil 
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rate of strain 

also frequency in the Fourier analysis 
also Stefan-Boltzmann's constant 

stress 

mean and fluctuating potential temperature 
soil logarithmic coordinates 

atmospheric log-linear coordinates 

zenith angle 

earth's angular velocity 

albedo 


also constant between $ and 1 used in the expression for 


the path length for water vapor 

also departure from isotropy in the ordering of terms 
various constants 

blackbody radiance 


heat capacities of dry air at constant pressure and at 
constant volume 


heat capacity of water vapor at constant pressure 


set of constants used in the computations of the surface 
temperature 


soil heat capacity 

coefficient of U(k,t +A t) and Vik yt +At) 
coefficient of G)(k,t,+A t) 

coefficient of Q(k,t,+A t) 

total derivative operator 

also half the soil depth in Myrup's (1969) model 


amplitude of the Fourier modes 
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derivative of & with respect to the soil linear coordinate 
second derivative of € with respect to Ze 

square of DDZ 

first derivative of Ywith respect to z 

ratio of the Y -derivative of f over 2 

second derivative of Y with respect to z 

square of DELZ 

constant used in various expressions 

first derivative of 2 with respect to ¥y 

first derivative of e ) with respect to Y 


diffusion operator 


eo 
also partial pressure of water vapor 


also coefficient in the resolvent cubic equation used to 
find the surface temperature 


latent heat flux 

also part of the solution of the quartic equation 

Coriolis parameter 

also used to represent any function 

also one of the coefficientsin the resolvent cubic equation 
computed value of a function F 

last iterative value of a function F 

new iterative value of F 

flux 


total infrared sky radiation received at the earth’s sur. 
face 


sensible heat flux from the atmosphore 
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J-i 

also an index 

net solar radiation 

intensity of radiation at frequency \) 

level at which the iterative formula computing the grid 
point height changes 

IKM + 4 

a general index 

wavenumber vector 

grid point 

von Karman'’s constant 

also a specific grid point 

coefficient of eddy diffusivity for momentum 
molecular soil thermal conductivity 

coefficient of eddy diffusivity for heat 

coefficient of eddy diffusivity for water vapor 
absorption coefficient of water vapor in a thin layer 


fourth-order tensor relating Os to the gradient of the 
mean wind J 


wavenumber 

also a general index 

mixing length 

maximum size of the eddies 

latent heat of vaporization 

number of increments in Taylor's series 


represents various exponents 
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top of the boundary layer in terms of k = fa) 
acceleration factor in the Gauss-Seidel iteration Lae, 
pressure 
reference pressure = 1000 mb 
turbulent Prandtl number 
fluctuating and mean mixing ratio for water vapor 
gas constant of the dry atmosphere 
gas constant of the water vapor 
infrared sky radiation 
terrestrial outgoing infrared radiation 
solar constant 
net radiation at the surface of the earth 
entropy 
entropy of dry air 
entropy of the water vapor 
soil heat flux 
time 
absolute temperature 
virtual temperature 
transmissivity of the cloudless atmosphere 
truncation error 
fluctuating and mean winds along the x direction 
friction velocity 
fluctuating and mean winds in the y direction 


geostrophic wind 
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w,W fluctuating and mean winds in the z direction 

Ww path length in water vapor 

x the first Cartesian coordinate in the horizontal 
y the second Cartesian coordinate in the horizontal 


also solution of the resolvent cubic equation 


Zz vertical Cartesian coordinate 
Z zenith solar angle 
Zo roughness height 
Zs proportion at each grid point of the total temperature 
difference between the surface and the bottom of the 
soil layer 


a roughness height 


24 canopy height in Myrup's (1969) model 


Zo top of the boundary layer in Myrup’s model 
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